CC BY
25
УДК 519.622 Е.А. Новиков
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ИОНИЗАЦИОННО-ДЕИОНИЗАЦИОННОГО ЦИКЛА ЦЕЗИЯ В ВЕРХНЕЙ АТМОСФЕРЕ L-УСТОЙЧИВЫМ МЕТОДОМ ВТОРОГО ПОРЯДКА
Описан алгоритм формирования правой части и матрицы Якоби дифференциальных уравнений химической кинетики. Численное моделирование цикла цезия в верхней атмосфере проведено посредством L-устойчивого метода второго порядка с контролем точности вычислений. Приведены результаты расчетов. Ключевые слова: цикл цезия, жесткая задача, 1-устойчивый метод, контроль точности.
Моделирование кинетики химических реакций применяется при исследовании разнообразных химических процессов. Предметом изучения являются временные зависимости концентраций реагентов, которые являются решением задачи Коши для системы обыкновенных дифференциальных уравнений. Трудности решения таких задач связаны с жесткостью и большой размерностью [1]. Здесь приведены результаты численного моделирования цикла цезия в верхней атмосфере Х-устойчивым методом второго порядка точности, в котором допускается замораживание как численной, так и аналитической матрицы Якоби [2].
Кинетическая схема химической реакции состоит из элементарных стадий, записываемых в виде [3]
N к N
Е “«Х -------’ (1)
1=1 1=1
где хI, 1</<Ыг - реагенты; &у, 1</<Ы - константы скоростей стадий; Ыг и N - соответственно, число реагентов и число стадий в реакции; агу и Ру, 1</<Ыт, 1</<Ы - стехиометрические коэффициенты. Процессу (1) в рамках сосредоточенной модели изотермического реактора постоянного объема соответствует система обыкновенных дифференциальных уравнений
С' = АТУ (2)
с заданным начальным условием С(0) = С0. Здесь АТ - стехиометрическая матрица, С и V - соответственно, вектор концентраций реагентов и скоростей стадий. В случае протекания реакции в не-
изотермических условиях к системе (2) добавляется уравнение теплового баланса
Т' = ^ТУ-а(Т - Т01)]/С^С, (3)
где Т - температура смеси в реакторе, Т01 - температура стенок реактора, О - вектор удельных теплот стадий, - вектор теплоемкостей реагентов, а = ая/г, а - коэффициент теплоотдачи, 5 и г -площадь поверхности и объем реактора соответственно. Верхний индекс Т у векторов О и СТр' означает транспонирование. Значения СТр' и а могут быть функциями концентраций реагентов, а а может еще зависеть от температуры.
Если реакция протекает в изотермическом реакторе постоянного объема с обменом вещества (открытая система, реактор идеального смешения), система дифференциальных уравнений записывается в виде
С' = АУТ + (Ср - С)/©, (4)
где Ср - вектор концентраций реагентов на входе, 0 - время пребывания смеси в реакторе, 0 = г/и, и - объемная скорость течения смеси через реактор. При протекании реакции в неизотермических условиях, система (4) дополняется уравнением теплового баланса
Т' = ^ТУ - а(Т - Т,і)]/С^С - (Т - Т02)/©, (5)
где Т02 - температура смеси на входе в реактор. Температура реагирующей смеси может задаваться в виде функции времени и концентраций.
Если стадия обратима, то за скорость стадии Ws принято считать разность скоростей прямого и обратного Ws- процессов, то есть Ws = Ws -Ws~, 1<5<^. Если в стадии участвует третья частица, то имеет место
N + N.
х=те, р = х ^,
і=1
где 8я, 1<7<^ - эффективности третьих частиц, ^п, 8я- и Сі, ^+1<1<^+^п - соответственно, количество, эффективности и концентрации инертных веществ. Значения компонент вектора Ws определяются из схемы химической реакции (1) по соотношениям
N + Міп ■
Ws+= к П О Ws-= к-, П с*, і=1 і=1
где к.. и к-., 1<.<К. - константы скоростей прямой и обратной стадий, соответственно. Они вычисляются по формуле к. = =AjTnJ ехр(-
Aj/RT), где Т - температура смеси в реакторе; Aj, п и Ej/R - заданные постоянные.
Стехиометрическая матрица АТ с элементами ау формируется из кинетической схемы (1) по следующему правилу. Номер стадии совпадает с номером столбца, а номер реагента с номером строки матрицы АТ. Если xi выступает как исходный реагент, то ау = -щ, если XI - продукт, то а^ = Ру. Если XI является одновременно исходным реагентом и продуктом, то ау = Ру-ау. Обычно в элементарной стадии участвует небольшое количество реагентов, то есть стехиометрическая матрица сильно разрежена.
Схема реакции ионизационно - деионизационного цикла цезия в верхней атмосфере имеет вид [4]: 1) 02'+С.+^С.+02, 2) С^+е^-С., 3) С. > С/+е, 4) 02+С+М^-СО2+М, 5) 02+е+М^0{+М, 6) 02 ^02+е, где к1 = 3.01010, к2=6.0-105, къ = 3.2410-3, к4 = 3.63 104, к5 = 3.63^104, к6 = 4.0^10-1. Реакция протекает с участием инертного вещества Ы2, причем концентрация N2] = 3.32-10"3. Эффективности третьих тел для всех реагентов равны 1, кроме эффективности 02 в пятой стадии, которая равна 12.4. Обозначим концентрации реагентов С1 = [е], С2 = [02-], С3 = [С.], С4 = [СО2], С5 = [С+], С6 = [О2]. Соответствующая система состоит из 6 дифференциальных уравнений вида
С1 = к2С1С5 + к3С3 - к5р2с1с6 + к6С2 , С2 = к1С2С5 + к5Р2С1С6 - к6С2 ,
С3 = к1С2С5 + к2С1С5 - к3С3 - к4Р1С3С6 , С4 = к4Р1С3С6 , (6)
С5 = -к1С2С5 - к2С1С5 + к3С3 , С6 = к1С2С5 - к4Р1С3С6 - к5Р2С1С6 + к6С2 ,
где ^1=С1+...+с6+^2], p2=Cl+_+c5+12.4c6+[N2]. Решение (6) осуществлялось на промежутке [0, 1000] с начальным шагом 10-5 при следующих начальных концентрациях реагентов: с1=1.66-10"16, С2=8.63-10"16, С3=1.6610-6, С4=0.0, С5=1.0310-15, С6=5.98^10-4. Сравнение эффективности алгоритма [2] проводилось с известным методом Гира dlsode [5] при точности расчетов е=10"2. Невысокая точность расчетов связана с тем, что метод [2] имеет второй порядок точности и поэтому проводить им высокоточные расчеты нецелесообразно. В качестве критерия эффективности выбраны if - число вычислений правой части и у - количество декомпозиций матрицы Якоби задачи (6) на интервале интегрирования. Алгоритму [2] для решения задачи (6) потребовалось 101 вычисление правой части и 14 декомпозиций матрицы Якоби. В методе dlsode требуемая точ-
ность 10-2 достигается при s=10-3 с вычислительными затратами //=192 и //=22.
Из результатов расчетов следует преимущество алгоритма [2] по числу вычислений правой части почти в два раза, в то время как преимущество по числу декомпозиций матрицы Якоби примерно в полтора раза. Наиболее эффективное применение данного алгоритма интегрирования предполагается на жестких задачах при точности расчетов s=10-2 (порядка 1%) и ниже.
Работа поддержана грантами РФФИ №08-01-00621 и Президента НШ-3431.2008.9
--------------------------------------------- СПИСОК ЛИТЕРАТУРЫ
1. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи / М.: Мир, 1999. - 685 с.
2. Новиков E.A. (2,1)-метод решения жестких неавтономных задач // Системы управления и информационные технологии. №2(32), 2008.-с. 12-15.
3. Эммануэль Н.М., Кнорре Д.Г. Курс химической кинетики / М.: Высшая школа, 1974. - 400с.
4. Edelson D. The new book in chemical kinetics // J. Chem. Ed., v. 52. 1975. -p.642.
5. http://www.netlib.org/odepackyindex.html nsrj=i
Novikov E.A
NUMERICAL SIMULATION IONIZATION-DEIONIZATION OF THE CYCLE OF CESIUM IN THE UPPER ATMOSPHERE OF L-STABLE METHOD OF SECOND ACCURACY ORDER
Algorithm right-hand side and Jacobian matrix of differential equations of chemical kinetics is described. Numerical simulation of the cycle of cesium in the upper atmosphere is realized by means of L-stable method of second accuracy order with the control accuracy. Results of computation are produced.
Key words: cycle of cesium, stiff problem, l-stability method, control accuracy.
— Коротко об авторе --------------------------------------------------
Новиков Евгений Александрович - доктор физико-математических наук, профессор, главный научный сотрудник Института вычислительного моделирования СО РАН (Красноярск),
E-mail: novikov@icm.krasn.ru
