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

V M. Vladimirov, E. F. Grinin, M. E. Sergiy, V N. Shepov

AUTOMATIC DEVICE FOR MEASURING RESISTIVITY OF SILICON FOUR-POINT PROBE METHOD

An automatic device for measuring the resistivity ofsingle-crystalline silicon by means of the four-point probe method had been developed.

Keywords: automatic device, single-crystalline silicon, resistivity, four-point probe method.

© Владимиров В. М., Гринин Э. Ф., Сергий М. Е., Шепов В. Н., 2009

УДК 519.622

А. Е. Новиков, Е. А. Новиков

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЦИКЛА ЦЕЗИЯ В ВЕРХНЕЙ АТМОСФЕРЕ ¿-УСТОЙЧИВЫМ МЕТОДОМ ВТОРОГО ПОРЯДКА ТОЧНОСТИ*

Описан алгоритм формирования правой части и матрицы Якоби дифференциальных уравнений химической кинетики. Численное моделирование цикла цезия в верхней атмосфере проведено посредством ¿-устойчивого метода второго порядка с контролем точности вычислений. Приведены результаты расчетов.

Ключевые слова: химическая кинетика, цикл цезия, жесткая задача, ¿-устойчивый метод, контроль точности.

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

Современные методы решения жестких задач используют обращение матрицы Якоби системы уравнений. В случае большой размерности исходной задачи декомпозиция такой матрицы практически полностью определяет общие вычислительные затраты. Для повышения эффективности расчетов в ряде алгоритмов используется замораживание матрицы Якоби, т. е. применение одной матрицы на нескольких шагах интегрирования [1;

2]. Наиболее успешно этот подход применяется в алгоритмах на основе многошаговых методов и, в частности, в формулах дифференцирования назад [3]. Однако использование данного подхода в алгоритмах интегрирования на основе известных безытерационных методов, к которым относятся методы типа Розенброка и их различные модификации, связано с определенными трудностями [1]. Ниже будет приведен алгоритм формирования правой части и матрицы Якоби дифференциальных уравнений химической кинетики, а также результаты численного моделирования ионизационно-де-ионизационного цикла цезия в верхней атмосфере ¿-устойчивым методом второго порядка точности, в ко-

тором допускается замораживание как численной, так и аналитической матрицы Якоби.

Кинетическая схема химической реакции состоит из элементарных стадий, записываемых в виде [4]

N к, ыг

'¿а,х1-------->^р,Х ’ (1)

/=1 1=1

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

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

С ’= ЛТУ (2)

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

T' = [QTV -a(T - Tm)]/CTVC,

(3)

где Т - температура смеси в реакторе; Т01 - температура стенок реактора; QT - вектор удельных теплот стадий; СТу - вектор теплоемкостей реагентов; а = аs/r, здесь а -коэффициент теплоотдачи; s и г - площадь поверхности и объем реактора соответственно. Верхний индекс Т у векторов QT и СТ¥ означает транспонирование. Теплоемкости реагентов и коэффициент теплоотдачи могут быть

* Работа поддержана грантами Российского фонда фундаментальных исследований № 08-01-00621 и Президента НШ-3431.2008.9

функциями концентраций реагентов, а коэффициент а может зависеть еще и от температуры.

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

Т дУ

Ь. = АТ —, 1 < ] < Ыг.

до,

(6)

С' = АКТ +1(Ср - С),

(4)

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

Т = ^ТУ - а(Т - Тт )]/СТС -©-(Т - Т>2), (5)

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

Если стадия обратима, то за скорость стадии Ж принимается разность скоростей прямого Ж+ и обратного Ж- процессов:

Ж = К + К, 1 < s < Ns.

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

ыг + ыт

V = Р-Ж, Р- = X е о, і =1

где е , 1 < і < N - эффективности третьих частиц; Ы.п, е и о, N + 1 < і < N + N - количество, эффективности и

гг г іп 7 т т

концентрации инертных веществ соответственно. Значения компонент вектора Ж определяются из схемы химической реакции (1) по соотношениям

Иг + Ит Ыг + Ит

К= к- П оГ-, К= к-* П оЬ% і=1 і=1

где к и к_, 1< - < N - константы скоростей прямой и обратной стадий соответственно, вычисляемые по формуле

кі = АТ"1 ехр(-Е. /ЯТ), здесь Т - температура смеси в реакторе, А ., п и Е/Я -заданные постоянные. Непосредственное использование этой формулы может приводить к неверному результату или переполнению арифметического устройства вследствие большого разброса данных постоянных [5; 6]. Поэтому для вычисления констант скоростей стадий используется следующее соотношение:

к] = ехр(1п А. + п] 1п Т - Е. /ЯТ).

Стехиометрическая матрица АТ с элементами а. формируется из кинетической схемы (1) по следующему правилу: номер стадии совпадает с номером столбца, а номер реагента - с номером строки матрицы АТ. Если х . выступает как исходный реагент, то а.. = -а.; если х . -продукт, то а . = р.; если же х . является одновременно исходным реагентом и продуктом, то а . = р.- а .

Обычно в элементарной стадии участвует небольшое количество реагентов, т. е. стехиометрическая матрица сильно разрежена. Тогда для системы (2) .-й столбец Ь. матрицы Якоби определяется по формуле

В случае совместного решения системы (2), (3) полученная матрица дополняется строкой ЬМг + 1 столбцом Ь ж и угловым элементом ЬМг+ 1^1г+1, которые определяются следующим образом:

Ь„, „..=- «?Т ^(Т - Т„,> +

г дс1 0С1

+[а(Т - Т01) - ^С / СТС} / СТС,

Ь,ж +1 = (Л , 1 <1< N.,

г дТ

т дV да т

Ь.г +., ,г +. = ^ — (Т - тт) - а] / СТУС -

. . дТ дТ

- ^ - а(Т - Т»)] (Щ;- С) / (СТ С )2. дТ

Для реактора с протоком диагональные элементы полученной матрицы подправляются с учетом добавки 1/0.

Нетрудно видеть, что компоненты вектора дV/дc. имеют вид

N. + N т

ду,/дС = а„Р,к,СТа-1 П С<л -

-Р^АС"-1 П сР,к - е, (ж; - ж;), (7)

к=1,к */

а для определения компонент вектора дV/дT используется следующие соотношения:

дк N' +

д», /дТ = Р, -Ь П <--дТ /=1

дк

••С-1=1 1 <' < *■ •

Если в ,-й стадии отсутствует третья частица, то в выражениях для ду /дс и ду /дТ нужно положить р = 1 и е = 0.

Использование (6) для определения элементов Ь,, 1 < 1 , , < N, матрицы Якоби позволяет применять формулу

Ь' = 1 • 1 < '•1 < •

Рассмотрим отдельное слагаемое (7), т. е.

а

р к а .оа-‘ 1 П оа-к -

-Г - - -і і X X і

-а . р к р 1 П с^л + а. е . (Ж + - Ж ).

к=1,к *1

Для определения данного выражения требуется выполнить три шага:

- на первом шаге формируется выражение а,,Р,к,а,1 и

а Р к р , при этом предполагается, что соотношения Р к

или Р к вычислены;

, —,

- на втором шаге определяются соотношения

са-- П

ста ' II с~“, ср 1 П ср-к'

- в случае присутствия в схеме реакции третьих частиц на третьем шаге задается выражение а,е,.(Ж,+-Ж,).

Опишем численный метод, который применялся для численного моделирования ионизационно-деионизаци-

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

у' = /('> у), у^ = у„, Ч < ' < 'к, (8)

где у и / - вещественные N-мерные векторы-функции; ' - независимая переменная, рассмотрим (т, к)-метод вида [7]

72

Уп+1 = Уп + акх + (1 - а)к2, а = 1 ——,

(9)

АА = М('п +р^ Уп )> АА = к1, где к1 и к2 - стадии метода; Бп = Е — акЛп; Е - единичная матрица; Ап - шаг интегрирования; Лп - некоторая матрица, представленная в виде Лп =/+ ЬБп + 0(к2); /п = д/у)ду -матрица Якоби системы (8); Бп - не зависящая от шага интегрирования матрица. Метод (9) можно применять с замораживанием как численной, так и аналитической матрицы Якоби.

Используя разложения стадий в ряды Тейлора, локальную ошибку 5п метода (9) запишем в виде

5п = (а -1/3)^ /;2 /+^ /;+1 /'у/1 +

+3 /;/ - 2 /у'/''-1 Б/ ] + 0(А4).

Тогда согласно [8] для контроля точности (9) используем оценку ошибки вида

еп = (а - 1)й2 // + 0(А3).

Учитывая, что

к2 - к1 = а^2/уд, + 0(А3Х

величину е, с точностью до членов 0(къ) оценим по формуле

е, = аТ1(а - 3)[к2 - ^].

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

V-, - • ^ ■■ - ае

е(.п ) = 11 А„-]" (к2 - к1) ||< 1/3

1/3 - а

где е - требуемая точность интегрирования; || • || - некоторая норма в К14; значение параметра,, 1 <, < 2, выбирается наименьшим, при котором выполняется данное неравенство. В конкретных расчетах норма 11 • || в этом неравенстве определяется по формуле

|[Б1;. (к2 - к)] 1|

Вп 1п (к2 -к1) || = шах-

где 1 - номер компоненты приближенного решения; 'г - положительный параметр. Если по 1-й компоненте решения выполняется неравенство |у' | < 'г, то контролируется абсолютная ошибка 'г • е, в противном случае -относительная ошибка е. В расчетах параметр 'г должен выбираться так, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой.

При численном вычислении матрицы Якоби шаг численного дифференцирования г., 1 <, < N, выбирается по формуле г. = тах(1014, 10-7|у.|). Решение жестких задач

обычно осуществляется с двойной точностью в силу плохой обусловленности матрицы Якоби системы (9). Константа 10-7 введена в формулу определения шага численного дифференцирования с целью его выдвижения на середину разрядной сетки. Теперь .-й столбец ап. матрицы Ап в формуле (9) вычисляется следующим образом:

/ (Ур..., У. + гУм ) - / (Ур..., У.Ум)

ап =------------------------------------------>

г

т. е. для определения матрицы Ап требуется N вычислений правой части задачи (8).

Попытка использования прежней матрицы Б осуществляется после каждого успешного шага интегрирования. Для того чтобы не ухудшить свойства устойчивости численной схемы, при замораживании матрицы Бп величина шага интегрирования тоже остается постоянной. Размораживание матрицы происходит в следующих случаях:

- когда нарушена точность расчетов;

- число шагов с замороженной матрицей достигло заданного максимального числа 1Н,

- прогнозируемый шаг больше последнего успешного в Qh раз;

- е(1) > е(2).

Числа 1к и Qh могут влиять на перераспределение вычислительных затрат. При 1к = 0 и Qh = 0 замораживания матрицы не происходит; при увеличении I и Qh число вычислений правой части задачи (8) возрастает, а количество обращений матрицы Якоби убывает. Эффективность алгоритма интегрирования также зависит от I и Q которые выбираются исходя из отношения стоимости вычисления функции / к стоимости декомпозиции матрицы Якоби и позволяют настраивать алгоритм на конкретный тип задач.

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

1) 02 + С- ® С- + 02

2) с;+ в ® с,

3) с- ® С+ + в

4) 02 + С + М ® С-02 + М

5) 02 + в + М ® 02- + М

6) 02 ® 02 + в

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

С1 = ^, С2 = [02 ^ С3 = [С-],

04 = [С-02], С5 = [С-+ ], 06 = [02].

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

с1 = к2С105 + к3с3 - к5Р2с!с6 + к6С2 ,

с2 = —k1c2c5 + к5p2c{6 — к6с2,

с’з = kiC2C5 + к2<сс5 — кзсз — к4PlC3C6 ,

с4 = кА Р1с3с6, (10)

с5 = —к1с2 с5 — к2с1с5 + к3с3,

с6 = к1с2с5 — к4 Р1сзсб — к5 Р2с1сб + кбс2.

где Р1 = с1 + с2 + с3 + с4 + с5 + с6 + [N2] ; р2 = с1 + с2 + с3 + + с4 + с5 + 12,4с6 +[N2] •

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

с1 = [e] = 1,66 -10—16, с2 = [O— ] = 8,63 -10—16,

с3 = [Cs ] = 1,66 -10—6,

с4 = CA] = 0,0, с5 = [C+ ] = 1,03 -10—15, с6 = [O2] = 5,98 -10—4.

Сравнение эффективности построенного алгоритма с известным методом Г ира dlsode в реализации А. Хинд-марша [10] при точности расчетов e = 10-2 показывает невысокую точность расчетов, которая связана с тем, что метод (9) имеет второй порядок точности и поэтому проводить с его помощью высокоточные расчеты нецелесообразно. В качестве критерия эффективности выбраны if - число вычислений правой части - и ij - количество декомпозиций матрицы Якоби задачи (10) на интервале интегрирования. Построенному алгоритму для решения задачи (10) потребовалось 101 вычисление правой части и 14 декомпозиций матрицы Якоби. В методе dlsode требуемая точность 10-2 достигается при e = 10-3 с вычислительными затратами if= 192 и ij = 22.

Таким образом, построенный алгоритм имеет преимущество по числу вычислений правой части почти в два раза и по числу декомпозиций матрицы Якоби - при-

мерно в полтора раза. Наиболее эффективное применение данного алгоритма интегрирования предполагается на жестких задачах при точности расчетов e = 10-2 (порядка 1 %) и ниже.

Библиографический список

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

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

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

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

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

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

7. Новиков, E. A. (2, 1)-метод решения жестких неавтономных задач / E. A. Новиков // Системы управления и информационные технологии. 2008. № 2 (32). С. 12-15.

8. Новиков, Е. А. Некоторые методы решения жестких систем, индуцированные одним и двумя вычислениями правой части / Е. А. Новиков, Ю. А. Шитов // Математические модели и методы решения задач механики сплошной среды. Красноярск, 1986. С. 11-19.

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

10. Brown, P. N. Reduced Storage Matrix Methods in Stiff ODE Systems / P. N. Brown, A. C. Hindmarsh // J. Appl. Math. & Comp. 1989. Vol. 31. P. 40-91.

A. E. Novikov, E. A. Novikov

NUMERUCAL SIMULATION OF THE CESIUM CYCLE IN THE UPPER ATMOSPHERE BY MEANS OF L-STABLE METHOD OF THE SECOND ORDER OF ACCURACY

An algorithm of right-hand side and Jacobian formation of differential equations of chemical kinetics is described. Numerical simulation of the cesium cycle in the upper atmosphere is conducted by means of the L-stable method of the second order of accuracy with the control accuracy. The results of the computation are presented.

Keywords: chemical kinetics, cycle of cesium, stiffproblem, L-stability method, control accuracy.

© Новиков А. Е., Новиков Е. А., 2009