close

Вход

Забыли?

вход по аккаунту

?

Применение метода Рунге-Кутта-Мерсона для решения диф ференциальных уравнений химической кинетики.

код для вставкиСкачать
Вестник КрасГАУ. 20 10. № 2
УДК 519.622
Л.В. Кнауб, A.E. Новиков, Ю.А. Шитов
ПРИМЕНЕНИЕ МЕТОДА РУНГЕ-КУТТА-МЕРСОНА ДЛЯ РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
ХИМИЧЕСКОЙ КИНЕТИКИ*
В статье описан алгоритм формирования правой части дифференциальных уравнений химической
кинетики. Численное моделирование цикла цезия в верхней атмосфере проведено посредством метода
Рунге-Кутта-Мерсона с контролем точности вычислений и устойчивости численной схемы.
Ключевые слова: химическая кинетика, цикл цезия, жесткая задача, методы Рунге-Кутта, контроль точности и устойчивости.
L.V. Knaub, A.Ye. Novikov, Yu.A. Shitov
THE RUNGE-KUTT-MERSON TECHNIQUE APPLICATION FOR SOLUTION OF THE CHEMICAL KINETICS
DIFFERENTIAL EQUATIONS
Algorithm for the formation of the right part of the chemical kinetics differential equations is described in the
article. Cesium cycle numerical modeling in the top atmosphere is made by means of the Runge-Kutt-Merson technique with the control of the calculations accuracy and numerical scheme stability.
Key words: chemical kinetics, cesium cycle, stiff problem, Runge-Kutt techniques, accuracy and stability
control.
Моделирование кинетики химических реакций применяется при исследовании разнообразных химических процессов [1]. Предметом изучения являются временные зависимости концентраций реагентов, которые составляют решение задачи Коши для системы обыкновенных дифференциальных уравнений. Трудности решения таких задач связаны с жесткостью и большой размерностью. Современные методы решения
жестких задач используют обращение матрицы Якоби системы уравнений [2–3]. В случае большой размерности исходной задачи, декомпозиция данной матрицы практически полностью определяет общие вычислительные затраты [4].
Здесь приведен алгоритм формирования правой части дифференциальных уравнений химической
кинетики. Даны результаты численного моделирования ионизационно-деионизационного цикла цезия в
верхней атмосфере явным методом Рунге-Кутта-Мерсона, в котором не используется обращение матрицы
Якоби. Наряду с точностью вычислений контролируется устойчивость численной схемы.
Кинетическая схема химической реакции состоит из элементарных стадий, записываемых в виде [1; 5–6]:
kj
Nr
Nr
x
x
ij i
ij i
i 1
(1)
i 1
где xi, 1≤i≤Nr – реагенты; kj, 1≤j≤Ns – константы скоростей стадий; Nr и Ns – соответственно число
реагентов и число стадий в реакции; αij и βij, 1≤i≤Nr, 1≤j≤Ns – стехиометрические коэффициенты.
Процессу (1) в рамках сосредоточенной модели изотермического реактора постоянного объема соответствует система обыкновенных дифференциальных уравнений
C
ATV
(2)
с заданным начальным условием C(0)=C0. Здесь AT – стехиометрическая матрица, C и V – соответственно
вектор концентраций реагентов и скоростей стадий. В случае протекания реакции в неизотермических условиях к системе (2) добавляется уравнение теплового баланса:
*Работа
поддержана грантами РФФИ №08-01-00621 и Президента НШ-3431.2008.9.
19
Математика и информатика
[QTV
T
(T T01 )] CVT C
(3)
где T – температура смеси в реакторе; T01 – температура стенок реактора; QT – вектор удельных теплот стадий; CTV – вектор теплоемкостей реагентов; α=άs/r, ά – коэффициент теплоотдачи; s и r – площадь
поверхности и объем реактора соответственно. Верхний индекс T у векторов QT и CTV означает транспонирование. Теплоемкости реагентов и коэффициент теплоотдачи могут быть функциями концентраций реагентов, а α может еще зависеть от температуры.
Если реакция протекает в изотермическом реакторе постоянного объема с обменом вещества (открытая система, реактор идеального смешения), система дифференциальных уравнений записывается в виде
1
AV T
C
C) ,
(C p
(4)
где Cp – вектор концентраций реагентов на входе; Θ – время пребывания смеси в реакторе; Θ=r/u, u –
объемная скорость течения смеси через реактор. При протекании реакции в неизотермических условиях
система (4) дополняется уравнением теплового баланса
T
[QTV
1
(T T01 )] CVT C
(T T02 ) ,
(5)
где T02 – температура смеси на входе в реактор. Температура реагирующей смеси может задаваться в
виде функции времени и концентраций, то есть T=T(t,C).
Если стадия обратима, то за скорость стадии Ws принято считать разность скоростей прямого Ws+ и
обратного Ws– процессов, то есть
Ws
Ws
Ns .
Ws , 1 s
Если в стадии участвует третья частица, то скорость перевычисляется по формуле
N r Nin
Vs
PW
s s
Ps
(6)
c
si i
i 1
где εsi, 1≤i≤Nr – эффективности третьих частиц, Nin, εsi и ci, Nr+1≤i≤Nr+Nin – соответственно количество
эффективности и концентрации инертных веществ. Значения компонент вектора Ws определяются из схемы
химической реакции (1) по соотношениям
N r Nin
Ws
ks
N r Nin
ci
si
Ws
i 1
k
ci si
s
i 1
где ks и k-s, 1≤s≤Ns – константы скоростей прямой и обратной стадий соответственно. Константы скоростей вычисляются по формуле
kj
n
AjT j exp( E j RT )
(7)
где T – температура смеси в реакторе; Aj, nj и Ej/R – заданные постоянные. Непосредственное использование данной формулы может приводить к неверному результату или переполнению арифметического
устройства вследствие большого разброса данных постоянных [5–6]. Поэтому для вычисления констант
скоростей стадий используется следующее соотношение:
20
Вестник КрасГАУ. 20 10. № 2
kj
exp(ln Aj
n j ln T
E j RT ) .
Стехиометрическая матрица AT с элементами αij формируется из кинетической схемы (1) по следующему правилу. Номер стадии совпадает с номером столбца, а номер реагента с номером строки матрицы AT.
Если xi выступает как исходный реагент, то aij=–αij, если xi – продукт, то aij=βij. Если xi является одновременно
исходным реагентом и продуктом, то aij=βij–αij. Обычно в элементарной стадии участвует небольшое количество реагентов, то есть стехиометрическая матрица сильно разрежена.
Опишем численный метод, который ниже будет применяться для численного моделирования ионизационно-деионизационного цикла цезия в верхней атмосфере. Для численного решения задачи Коши для
системы обыкновенных дифференциальных уравнений
y
f (t y)
y(t0 )
y0
(8)
t0 t tk
где y и f – вещественные N-мерные вектор-функции; t – независимая переменная. Рассмотрим метод
Рунге-Кутта-Мерсона [7]:
1
2
1
k
k
k5
1
1
4
6
3
6
1
k1 ) k3 hf ( yn
k1 hf ( yn ) k2 hf ( yn
3
1
1
3
k1
k4 hf ( yn
k1
k3 ) k5 hf ( yn
8
8
2
yn
yn
1
k1
6
3
k3
2
1
k2 )
6
(9)
2k 4 )
Пятое вычисление правой части не дает увеличение порядка до пятого, но позволяет расширить интервал устойчивости до 3,5, и оценить локальную ошибку δn,4 с помощью ранее вычисленных стадий ki,
1≤i≤5, то есть
n4
(2k1 9k3 8k4
k5 ) / 30
Теперь можно записать неравенство для контроля точности через оценку локальной ошибки, то есть
n4
5
5 4
(10)
Одна из причин успеха схемы (9) заключается по-видимому в экономичном способе оценки ошибки, на
основе которой осуществляется контроль точности вычислений и автоматический выбор величины шага интегрирования. Кроме того, область устойчивости метода (9) достаточно велика как по вещественной, так и по
мнимой оси комплексной плоскости. Это позволяет использовать его при решении жестких задач, собственные числа матрицы Якоби которых имеют достаточно большую мнимую часть.
Как показывают расчеты, контроль устойчивости численной схемы приводит не только к повышению
эффективности алгоритма интегрирования, но и позволяет вычислить решение ряда задач, которые не удается рассчитать алгоритмом без контроля устойчивости. Поэтому построим неравенство для контроля устойчивости метода (9). Применяя к разности (k3–k2) формулу Тейлора с остаточным членом в форме Лагранжа первого порядка, будем иметь
k3 k2
h f ( y n)
( k2
6
y
21
k1 )
Математика и информатика
где вектор
y n
вычислен в некоторой окрестности решения
k2
k1
(h 2 / 3) f n f n
y(tn ) . Учитывая, что
O(h3 )
для контроля устойчивости (9) можно использовать неравенство
6max | (k3 k2 )i /(k2 k1 )i | 3 5
(11)
1 i N
где числу 3,5 примерно равна длина интервала устойчивости численной схемы. Применяя метод (9)
для решения задачи
y
y0 , Re( ) 0 ,
y, y(0)
получим
yn
1
Q( x) yn , x h ,
Q( x ) 1 x
1 2
x
2
1 3
x
6
1 4
1 5
x
x.
24
120
Область устойчивости метода приведена на рисунке. Отметим, что по вещественной и мнимой оси
область устойчивости ограничена числом 3,5.
Область устойчивости метода (9)
Сформулируем алгоритм переменного шага с контролем точности вычислений и устойчивости численной схемы. Пусть приближение к решению yn в точке tn вычислено с шагом интегрирования hn .
1. Вычисляется стадия k1 по формуле (9).
2. Вычисляются стадии k2, k3, k4 и k5 по формуле (9).
3. Вычисляется оценка ошибки εn,4 по формуле εn,4=||2k1–9k3+8k4–k5||/150.
4. Вычисляется число q1 по формуле q15εn,4=ε, где ε есть требуемая точность интегрирования. При определении q1 учитывается, что εn,4=O(hn5).
5. Если q1<1, то есть требуемая точность не выполняется, то hn полагается равным q1hn и происходит возврат на шаг 2.
22
Вестник КрасГАУ. 20 10. № 2
6. Вычисляется приближенное решение yn+1 по формуле (9).
7. Вычисляется значение tn+1=tn+hn.
8. Вычисляется оценка максимального собственного числа vn,4 матрицы Якоби по формуле
vn,4=6max1≤i≤N|k3i–k2i|/|k2i–k1i|.
9. Вычисляется число q2 по формуле q2vn,4=3,5. При определении q2 учитывается, что vn,4=O(hn). Отметим, что число q2 ограничивает размер шага по устойчивости.
10. Если q2<q1, то есть шаг по устойчивости меньше шага по точности, то новый шаг интегрирования
hn+1 выбирается равным hn. В противном случае, hn+1 вычисляется по формуле =min(q1, q2)hn.
11. Выполняется следующий шаг интегрирования по методу четвертого порядка точности, то есть
происходит переход на шаг 1.
Ниже алгоритм интегрирования на основе численной схемы (9), в котором для контроля точности и
устойчивости используются соответственно неравенства (10) и (11) будем называть RK4ST.
Теперь рассмотрим модель ионизационно-деионизационного цикла цезия в верхней атмосфере. Настоящая модель извлечена из большой кинетической схемы и широко используется за рубежом как типичный пример жестких систем кинетики [8]. Схема реакции имеет вид:
1)
O2
Cs
Cs
3)
Cs
Cs
e,
5)
O2
e M
O2 ,
4) O2
O2
2)
Cs
Cs ;
e
CsO2 M ;
O2 e ,
Cs M
M , 6) O2
где константы скоростей стадий имеют вид: k1=3,0∙1010, k2=6,0∙105, k3=3,24∙10-3, k4=3,63∙104,
k5=3,63∙104, k6=4,0∙10-1. Реакция протекает с участием инертного вещества N2, причем концентрация
[N2]=3,32∙10-3. Эффективности третьих тел для всех реагентов равны 1, кроме эффективности O2 в пятой
стадии, которая равна 12,4. Обозначим концентрации реагентов следующим образом:
c1 [e] , c2 [O2 ] , c3 [Cs ] ,
c4 [CsO2 ] , c5 [Cs ] , c6 [O2 ] .
Соответствующая система состоит из 6 дифференциальных уравнений вида:
c1
k2c1c5 k3c3 k5 p2c1c6 k6c2 ,
c2
k1c2c5 k5 p2c1c6 k6c2 ,
c3 k1c2c5 k2c1c5 k3c3 k4 p1c3c6 ,
c4 k4 p1c3c6 ,
c5
k1c2c5 k2c1c5 k3c3 ,
c6 k1c2c5 k4 p1c3c6 k5 p2c1c6 k6c2 ,
(12)
где
p1 c1 c2 c3 c4 c5 c6 [ N2 ] ,
p2 c1 c2 c3 c4 c5 12.4c6 [ N2 ] .
Интегрирование осуществлялось на промежутке [0,1000] с начальным шагом 10-5 при следующих начальных концентрациях реагентов:
c1 [e] 1 66 10
16
,
c2
23
[O2 ] 8 63 10
16
,
Математика и информатика
c3 [Cs ] 1 66 10
c5 [Cs ] 1 03 10
15
6
,
,
c4 [CsO2 ] 0 0 ,
c6 [O2 ] 5 98 10 4 .
Сравнение эффективности построенного алгоритма проводилось с методом Мерсона (MERSON) без
контроля устойчивости при точности расчетов ε=10-6. В качестве критерия эффективности выбрано число
вычислений правой части задачи (12) на интервале интегрирования. Алгоритму RK4ST с контролем точности
вычислений и устойчивости численной схемы для решения задачи (12) потребовалось примерно в полтора
раза меньше вычислений правой части по сравнению с алгоритмом MERSON. Это является следствием контроля устойчивости численной схемы, который позволяет устранить неоправданные возвраты. Отметим, что
при расчетах жестких задач явными методами без контроля устойчивости почти каждый шаг сопровождается
повторным вычислением решения из-за возникающей неустойчивости, что приводит к понижению эффективности расчетов.
Наиболее эффективное применение данного алгоритма интегрирования предполагается на нежестких
задачах, а также задачах средней жесткости при точности расчетов ε порядка 10-4~10-6.
Литература
1.
2.
3.
4.
5.
6.
7.
8.
Кнорре Д.Г., Эммануэль Н.М. Курс химической кинетики. – М.: Высш. шк., 1974. – 400 с.
Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. – М.: Мир, 1999. – 685 c.
Новиков Е.А. Явные методы для жестких систем. – Новосибирск: Наука, 1997. – 197с.
Byrne G.D., Hindmarsh A.C. ODE solvers: a review of current and coming attractions // J. of Comput. Physics. – 1987. – №70. – P. 1 – 62.
Babusok V.I., Novikov E.A. Numerical solution of direct kinetic problems // React. Kinet. Catal. Lett. – 1982. –
Vol. 21. – P.121–124.
Новиков Е.А., Бабушок В.И. Генератор правой части и матрицы Якоби дифференциальных уравнений
химической кинетики. – Новосибирск: ВЦ СО АН СССР, 1982. – 27с.
Merson R.H. An operational methods for integration processes // Proc. of Symp. on Data Processing. Weapons Research Establishment. – Salisbury, Australia, 1957. – P.239 – 240.
Edelson D. The new book in chemical kinetics // J. Chem. Ed. – 1975. – Vol. 52. – P. 642–643.
24
Вестник КрасГАУ. 20 10. № 2
25
Документ
Категория
Без категории
Просмотров
8
Размер файла
620 Кб
Теги
рунге, мерсона, решение, уравнения, кутта, ференциальных, метод, кинетике, диф, применению, химические
1/--страниц
Пожаловаться на содержимое документа