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

УДК 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]:

ЛГ. к] ТУ,.

2Х*#------------->2^РуХп (1)

7=1 7=1

где Xi, 1<i<Nr - реагенты; к/, 1</<Ns - константы скоростей стадий; Nr и Ns - соответственно число

реагентов и число стадий в реакции; а/ и р/, 1</<М-, 1</<Ns - стехиометрические коэффициенты.

Процессу (1) в рамках сосредоточенной модели изотермического реактора постоянного объема соответствует система обыкновенных дифференциальных уравнений

С = АТУ (2)

с заданным начальным условием С(0)=Со. Здесь Ат- стехиометрическая матрица, С и V - соответственно вектор концентраций реагентов и скоростей стадий. В случае протекания реакции в неизотермических условиях к системе (2) добавляется уравнение теплового баланса:

*Работа поддержана грантами РФФИ №08-01-00621 и Президента НШ-3431.2008.9.

Т' = \<2ТУ - а(Т - Т01)\/СуС, (3)

где Т - температура смеси в реакторе; Т01 - температура стенок реактора; 0Т - вектор удельных теп-лот стадий; 0Ту - вектор теплоемкостей реагентов; а=ав1г, а - коэффициент теплоотдачи; 5 и г - площадь поверхности и объем реактора соответственно. Верхний индекс Т у векторов 0Т и 0Ту означает транспонирование. Теплоемкости реагентов и коэффициент теплоотдачи могут быть функциями концентраций реагентов, а а может еще зависеть от температуры.

Если реакция протекает в изотермическом реакторе постоянного объема с обменом вещества (открытая система, реактор идеального смешения), система дифференциальных уравнений записывается в виде

г 1

С = А¥т+-(Ср-С), (4)

(у

где 0Р - вектор концентраций реагентов на входе; 0 - время пребывания смеси в реакторе; 0=г/и, и -объемная скорость течения смеси через реактор. При протекании реакции в неизотермических условиях система (4) дополняется уравнением теплового баланса

Г = [ОтГ-а(Т-Тт )]/СуС -1(т-т02), (5)

0

где Т02 - температура смеси на входе в реактор. Температура реагирующей смеси может задаваться в виде функции времени и концентраций, то есть Т=Т(?,0).

Если стадия обратима, то за скорость стадии Ws принято считать разность скоростей прямого и обратного Ws~ процессов, то есть

Если в стадии участвует третья частица, то скорость перевычисляется по формуле

мг+м1п

К=РК, Р,= Це„с„ (6)

1=1

где £5/, 1</<Ы - эффективности третьих частиц, Ы/п, г» и с/, Ыг+1</<Ыг+Ып - соответственно количество эффективности и концентрации инертных веществ. Значения компонент вектора Ws определяются из схемы химической реакции (1) по соотношениям

К=К П с°",ж; = к_, П с?-,

2=1 2=1

где 1<5 и к.*, 1<5<Ы - константы скоростей прямой и обратной стадий соответственно. Константы скоростей вычисляются по формуле

к} = Ар ехр(-Я/КГ), (7)

где Т - температура смеси в реакторе; Д, п и Е^ - заданные постоянные. Непосредственное использование данной формулы может приводить к неверному результату или переполнению арифметического

устройства вследствие большого разброса данных постоянных [5-6]. Поэтому для вычисления констант скоростей стадий используется следующее соотношение:

kJ - exp(ln Aj + rij InT - Ej/RT).

Стехиометрическая матрица AT с элементами a,j формируется из кинетической схемы (1) по следующему правилу. Номер стадии совпадает с номером столбца, а номер реагента с номером строки матрицы AT. Если X, выступает как исходный реагент, то aj=-aj, если x/ - продукт, то aj=pj. Если x/ является одновременно исходным реагентом и продуктом, то Sii=eij-aij. Обычно в элементарной стадии участвует небольшое количество реагентов, то есть стехиометрическая матрица сильно разрежена.

Опишем численный метод, который ниже будет применяться для численного моделирования иониза-ционно-деионизационного цикла цезия в верхней атмосфере. Для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений

y' = f(t,y\ У(Ч) = у0, tQ<t<tk, (8)

где у и f- вещественные W-мерные вектор-функции; t - независимая переменная. Рассмотрим метод Рунге-Кутта-Мерсона [7]:

!/ 2i li

Уп+1=Уп+Тк1+тК+тк5,

6 3 6

К = ¥(УпХ К = ¥(у„ + \ю, къ = hf (уп + hix + -U2), (9)

3 6 6

К = + к5 =hf(y„ +1*.-h, + 2к,).

Пятое вычисление правой части не дает увеличение порядка до пятого, но позволяет расширить интервал устойчивости до 3,5, и оценить локальную ошибку 6„,4 с помощью ранее вычисленных стадий к/, 1</<5, то есть

Sn4= (2кх -9к3 + 8к4 -к5)/30.

Теперь можно записать неравенство для контроля точности через оценку локальной ошибки, то есть

|£'й,4|<5^/4. (10)

Одна из причин успеха схемы (9) заключается по-видимому в экономичном способе оценки ошибки, на

основе которой осуществляется контроль точности вычислений и автоматический выбор величины шага интегрирования. Кроме того, область устойчивости метода (9) достаточно велика как по вещественной, так и по мнимой оси комплексной плоскости. Это позволяет использовать его при решении жестких задач, собственные числа матрицы Якоби которых имеют достаточно большую мнимую часть.

Как показывают расчеты, контроль устойчивости численной схемы приводит не только к повышению эффективности алгоритма интегрирования, но и позволяет вычислить решение ряда задач, которые не удается рассчитать алгоритмом без контроля устойчивости. Поэтому построим неравенство для контроля устойчивости метода (9). Применяя к разности (ta-fc) формулу Тейлора с остаточным членом в форме Лагранжа первого порядка, будем иметь

, h д/( у )

3 2 = ^ ~ ( 2 1)’

6 ду

где вектор у вычислен в некоторой окрестности решения у{^п). Учитывая, что

к2-к1=(ь2/з)г/п+о(ь3),

для контроля устойчивости (9) можно использовать неравенство

бтах | (к3 - к2)і /{к2 - к1)і |< 3.5,

(11)

где числу 3,5 примерно равна длина интервала устойчивости численной схемы. Применяя метод (9) для решения задачи

Область устойчивости метода приведена на рисунке. Отметим, что по вещественной и мнимой оси область устойчивости ограничена числом 3,5.

Сформулируем алгоритм переменного шага с контролем точности вычислений и устойчивости численной схемы. Пусть приближение к решению уп в точке ^ вычислено с шагом интегрирования Ип.

1. Вычисляется стадия к1 по формуле (9).

2. Вычисляются стадии кг, кз, к4 и кь по формуле (9).

3. Вычисляется оценка ошибки спа по формуле £„,4=||2к1-9кз+8к4-к5||/150.

4. Вычисляется число Ц1 по формуле д15£п,4=£, где е есть требуемая точность интегрирования. При определении Ц1 учитывается, что £„,4=0(Л„5).

5. Если ф<1, то есть требуемая точность не выполняется, то ^ полагается равным д1Лп и происходит возврат на шаг 2.

у' = Лу, у(0) = уо, Яе(Я)<0,

получим

Область устойчивости метода (9)

6. Вычисляется приближенное решение уп+1 по формуле (9).

7. Вычисляется значение и+1=^„.

8. Вычисляется оценка максимального собственного числа Vn,4 матрицы Якоби по формуле

^,4=6таХ 1< йу| к3к2 '|/| к2к1'|.

9. Вычисляется число Ц2 по формуле q2v„4=3,5. При определении q2 учитывается, что уп,4=0(Ь„). Отметим, что число q2 ограничивает размер шага по устойчивости.

10. Если q2<q1, то есть шаг по устойчивости меньше шага по точности, то новый шаг интегрирования Л„+1 выбирается равным Ь„. В противном случае, Л„+1 вычисляется по формуле =т1п(ф, q2)hn.

11. Выполняется следующий шаг интегрирования по методу четвертого порядка точности, то есть происходит переход на шаг 1.

Ниже алгоритм интегрирования на основе численной схемы (9), в котором для контроля точности и устойчивости используются соответственно неравенства (10) и (11) будем называть РК4ЭТ.

Теперь рассмотрим модель ионизационно-деионизационного цикла цезия в верхней атмосфере. Настоящая модель извлечена из большой кинетической схемы и широко используется за рубежом как типичный пример жестких систем кинетики [8]. Схема реакции имеет вид:

1 )0~2+с;^ся+02,2) с; I ; з) с, —>■ с; + е, 4)02+Сл, + м —> ср2 + м ]

5) 02 + 6 + А/ —^ 02 + А/, 6) 02 —^ 02 + 6,

где константы скоростей стадий имеют вид: к1=3,01010, к2=6Д105, кз=3,24^10-3, к4=3,63^104, к5=3,63^104, кб=4,0^10-1. Реакция протекает с участием инертного вещества N2, причем концентрация [^2]=3,32-10-3. Эффективности третьих тел для всех реагентов равны 1, кроме эффективности О2 в пятой стадии, которая равна 12,4. Обозначим концентрации реагентов следующим образом:

сх — \в\, с2 — \02 ], с3 — [С5 ],

с4=[с,о2], с5 = [с;], с6=[о2].

Соответствующая система состоит из 6 дифференциальных уравнений вида:

С\ — ~к2С\СЪ кзсз ~ к5Р2^\^6 к6С2 ’

с2 — —кхс2с5 + к5р2схс6 — к( с2,

с3 = к]с2с5 + к2схс5 — к3с3 - кАрхс3сь,

с\ = кАрхс3сь, (12)

с'5 = -кхс2с5 — к2схс5 + к3с3,

с6 — кхс2с5 — кАрхс3сь — к5р2схс6 + к(с2,

где

рх=с х+с2+с3+сА+с5+с6 + [А^2], р2 = сх + с2 + с3 + сА + + 12.4с6 + [Л^] ■

Интегрирование осуществлялось на промежутке [0,1000] с начальным шагом 10-5 при следующих начальных концентрациях реагентов:

Сх = [е\ = 1.66 • 10 16, с2 = [02] = 8.63 • 10 16,

съ = [CJ = 1.66-10б, с4 = [Ср2] = 0.0, с5 = [с;] = 1.03 10 15, с6 =[02] = 5.98• 10 4.

Сравнение эффективности построенного алгоритма проводилось с методом Мерсона (MERSON) без контроля устойчивости при точности расчетов £=10"6. В качестве критерия эффективности выбрано число вычислений правой части задачи (12) на интервале интегрирования. Алгоритму RK4ST с контролем точности вычислений и устойчивости численной схемы для решения задачи (12) потребовалось примерно в полтора раза меньше вычислений правой части по сравнению с алгоритмом MERSON. Это является следствием контроля устойчивости численной схемы, который позволяет устранить неоправданные возвраты. Отметим, что при расчетах жестких задач явными методами без контроля устойчивости почти каждый шаг сопровождается повторным вычислением решения из-за возникающей неустойчивости, что приводит к понижению эффективности расчетов.

Наиболее эффективное применение данного алгоритма интегрирования предполагается на нежестких задачах, а также задачах средней жесткости при точности расчетов £ порядка 10-4~10-6.

Литература

1. Кнорре Д.Г., Эммануэль Н.М. Курс химической кинетики. - М.: Высш. шк., 1974. - 400 с.

2. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. - М.: Мир, 1999. - 685 с.

3. Новиков Е.А. Явные методы для жестких систем. - Новосибирск: Наука, 1997. - 197с.

4. Byrne G.D., Hindmarsh A.C. ODE solvers: a review of current and coming attractions // J. of Comput. Physics. - 1987. - №70. - P. 1 - 62.

5. Babusok V.I., Novikov E.A. Numerical solution of direct kinetic problems // React. Kinet. Catal. Lett. - 1982. -Vol. 21. - P. 121-124.

6. Новиков Е.А., Бабушок В.И. Генератор правой части и матрицы Якоби дифференциальных уравнений химической кинетики. - Новосибирск: ВЦ СО АН СССР, 1982. - 27с.

7. Merson R.H. An operational methods for integration processes // Proc. of Symp. on Data Processing. Weapons Research Establishment. - Salisbury, Australia, 1957. - P.239 - 240.

8. Edelson D. The new book in chemical kinetics // J. Chem. Ed. - 1975. - Vol. 52. - P. 642-643.

'--------♦------------