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

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

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

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

N r ^k j N r

5 a ijxi—>5irx-          (1)

i =1                              i =1

где x_i , 1 ≤ i ≤ N_r – реагенты; k_j , 1 ≤ j ≤ N_s – константы скоростей стадий; N_r и N_s – число реагентов и число стадий в реакции соответственно; a _ij и β _ij , 1 ≤ i ≤ N_r , 1 ≤ j ≤ N_s – стехиометрические коэффициенты.

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

C ′= A^TV                   (2)

с заданным начальным условием C (0) = C ₀. Здесь A^T – стехиометрическая матрица, C и V – векторы концентраций реагентов и скоростей стадий соответственно. В случае протекания реакции в неизотермических условиях к системе (2) добавляется уравнение теплового баланса

T ′= [ Q^TV - a ( T - T 01 )] / C V ^TC ,             (3)

где T – температура смеси в реакторе; T01 – температура стенок реактора; QT – вектор удельных теплот стадий; CTV – вектор теплоемкостей реагентов; a = as/r, здесь a – коэффициент теплоотдачи; s и r – площадь поверхности и объем реактора соответственно. Верхний индекс T у векторов QT и CTV означает транспонирование. Теплоемкости реагентов и коэффициент теплоотдачи могут быть функциями концентраций реагентов, а коэффициент a может зависеть еще и от температуры.

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

C ‘ = AV^T + !( C p - C ) ,              (4)

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

T ‘ = [ Q T V - a ( T - T ₀₁)] Z C C - 1 ( T - T ₀ 2 ) ,     ⁽⁵⁾

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

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

W s = W s+ + W s - , 1 <  “ <  N s .

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

N r + N in

У =PW, Р = У £ с., sss,ssii, i=1

где s si , 1 <  i <  N r - эффективности третьих частиц; N in , £ si и c i , N r + 1 <  i <  N r + N in - количество, эффективности и концентрации инертных веществ соответственно. Значения компонент вектора W_s определяются из схемы химической реакции (1) по соотношениям

N r + N in                       N r + N in

W+=kcasi,W-=kcβsi, ssis-si i =1                                     i =1

где к и k , 1< s < N - константы скоростей прямой и обратной стадий соответственно, вычисляемые по формуле kj = AjTn exp(-EjIRT), здесь T - температура смеси в реакторе, Aj, nj и Ej/R -заданные постоянные. Непосредственное использование этой формулы может приводить к неверному результату или переполнению арифметического устройства вследствие большого разброса данных постоянных [5; 6]. Поэтому для вычисления констант скоростей стадий используется следующее соотношение:

k j = exp(ln A j + n j ln T - E jI RT ).

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

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

В случае совместного решения системы (2), (3) полученная матрица дополняется строкой b Nr _{+ 1} . , столбцом b i. _{Nr +1} и угловым элементом b N, ₊ 1, N , ₊₁, которые определяются следующим образом:

bN+1,i= Q dV-|a (T - T01) + r            ∂ci ∂ci

+ [ a ( T - T 01 ) - QV ] cj C V C } Z C V C ,

ЗУ b., .,+1= (A dT).,

r ^,

b_{N+1 w+}, = [ QT — -— ( T - T ,) -a ] I CTC - N r + 1 , N r + 1 L     d T d T         01 -I V

T

- [ QV - a ( T - T )](— C ) I ( C V C ) ².

∂ T

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

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

Nr+Nin д vs Id ci = asipskscc' “ 1 П c. “ -k=1, k * i

N_r + N_in

-e “ P s k - s c “ ^{- 1} П c i “k - £ si (W “ + - W “ ^- ) ,        (7)

к=1, k * i а для определения компонент вектора д V/дTиспользует-ся следующие соотношения:

N r + N in

∂v/∂T=p∂ksc

“          д T   =1

_{∂ k} N R + N in

- P_s —- П c в “ , 1 <  “ <  N .

“ д t 1=1 i“

Если в s -й стадии отсутствует третья частица, то в выражениях для д v_s / д c i и д v_s / д T нужно положить p_s = 1 и s _si = 0.

Использование (6) для определения элементов bj, 1 < i, j < Nr, матрицы Якоби позволяет применять формулу b Za- fv^, 1 < i,j < N,.

s =1       ^д c i

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

Nr+Nin ajsPsks asica “i-1 П c.a “k - k=1, k * i

N r + N in

-a,$p$k $в$,cf“i-1 П cв“k + a,£si(W++-WГ)• jss-ssiiijssiss k=1, k * i

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

• -    на первом шаге формируется выражение a_jsp_sk_s a _si и a j P_sk_-s в _si , при этом предполагается, что соотношения p_sk_s или pk_s вычислены;

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

  N r + N in

  c i " ^{- 1} П c i

  k =1, k * i

  
  

Nr+Nin c. s.-1 П  ci “k;

k =1, k * i

• -    в случае присутствия в схеме реакции третьих частиц на третьем шаге задается выражение a_js s _si ( W_s ⁺- W_s ^-).

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

У ' = f ^{( t} , У ) , У ^{( t} 0 ) = У 0 , ^t 0 <  ^t <  t k ,           (8)

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

, „ „        , V 2

Уп+1 = Уп + ak1 + (1 - a)к2 , a = 12",

Dnki = hnf (tn +в h, Уп), Dnk 2 = ki, где к 1 и к2 - стадии метода; Dn = E - ahA; E - единичная матрица; hn - шаг интегрирования; An - некоторая матрица, представленная в виде An =fп+ hBn + O(h2); fn = dfyn)/5у -матрица Якоби системы (8); Bn - не зависящая от шага интегрирования матрица. Метод (9) можно применять с замораживанием как численной, так и аналитической матрицы Якоби.

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

5 п = ( a - 1 / 3) h³fy'² f + h/ '+ ¹ f y f ² +

24     6

+ 3 f t f - 2 f ; f, '- 2 B n f ] + о ( h ⁴) .

Тогда согласно [8] для контроля точности (9) используем оценку ошибки вида sn = (a -3)h2 f‘f + O(h3).

Учитывая, что к2 - ki = ah2fy',nfn + O(h3), величину sn с точностью до членов O(h3) оценим по формуле sп = a 1(a - з)[к2 - ki].

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

s ( У п ) = I| D п" j" ( к 2 - к Д||< —— ,

1/3 - a где s - требуемая точность интегрирования; || ■ || - некоторая норма в RN; значение параметра у", 1 < у" < 2, выбирается наименьшим, при котором выполняется данное неравенство. В конкретных расчетах норма || ■ || в этом неравенстве определяется по формуле

II D    ( к 2

-

к ,) || = max

¹          1< i <  N

| [ D^j (к2 - кды I у" | + tr где i - номер компоненты приближенного решения; tr - положительный параметр. Если по i-й компоненте решения выполняется неравенство yi | < tr, то контролируется абсолютная ошибка tr ■ s, в противном случае -относительная ошибка s. В расчетах параметр tr должен выбираться так, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой.

При численном вычислении матрицы Якоби шаг численного дифференцирования rj, 1 < j < N, выбирается по формуле rj = max(10-14, 10 7y|). Решение жестких задач обычно осуществляется с двойной точностью в силу плохой обусловленности матрицы Якоби системы (9). Константа 10-7 введена в формулу определения шага численного дифференцирования с целью его выдвижения на середину разрядной сетки. Теперь j-й столбец a" матрицы A" в формуле (9) вычисляется следующим образом:

a

п =

^{f (} У 1 ,-, У у ⁺ r j,-, У м ) ^{- f (} У 1

., У у

У м )

r j

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

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

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

• -    число шагов с замороженной матрицей достигло заданного максимального числа I h ;

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

• -    s (1) >  s (2).

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

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

• 1)    о - + с ; >  c s + o ₂

• 2)    C ;+ e ^ C_s

• 3)    C_s ^ с ;+ e

• 4)    о ₂ + C_s + M ^ C_sO₂ + M

• 5)    O ₂ + e + M ^ O ₂ ^- + M

• 6)    о ₂ ^- ^ O ₂ + e

  где константы скоростей стадий имеют вид к ₁ = 3,0 ■ 10¹⁰, к ₂= 6,0 ■Ю⁵, к ₃= 3,24 ■ 10^-3, к ₄= 3,63 ■ 10⁴, к ₅= 3,63 ■ 10⁴, к ₆= 4,0 ■ 10^-1. Реакция протекает с участием инертного вещества N ₂, причем концентрация [ N ₂] = 3,32 ■ 10^-3. Эффективности третьих тел для всех реагентов равны 1, кроме эффективности о ₂ в пятой стадии, которая равна 12,4. Обозначим концентрации реагентов следующим образом:

С1 = [e], c2 = [O;], c3 = [Cs], c4 = [CsO2], c5 = [cs], c6 = [O2].

Соответствующая система состоит из шести дифференциальных уравнений вида c1 = - к2 c1 c 5 + к3 c3 - к5 p 2 c1 c 6 + к6 c 2, c 2 = - k1 c 2 c5 + k5 P 2 C1 c 6 - k6 c 2 , cз = ki c2 c5 + k2 ci c5 — kз cз — k4 Pi cз c6, c4= k4 Pic3 c 6,                   (10)

c 5 = — ki c 2 c5 — k 2 ci c5 + k3 c 3 , c 6 = ki c 2 c 5 — k4 pi c3 c 6 — k5 p 2 ci c6 + k6 c 2’ где pi = ci + c2 + c3 + c4 + c5 + c6 + [ N2 ]; p2 = ci + c2 + c3 + + c4 + c5 +i2,4c6 + [N2].

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

c _i = [ e ] = i,66 - i0 ^{— i6}, c 2 = [ O ₂ ^— ] = 8,63 - i0^—16, c 3 = [ C s ] = i,66 - I0 ^{— 6}, c 4 = [ C s O 2 ] = 0,0, c 5 = [ C s ] = i,03 - i0 ^{— i5}, c 6 = [ O ₂] = 5,98 - I0 ^- 4 .

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

Таким образом, построенный алгоритм имеет преимущество по числу вычислений правой части почти в два раза и по числу декомпозиций матрицы Якоби – при- мерно в полтора раза. Наиболее эффективное применение данного алгоритма интегрирования предполагается на жестких задачах при точности расчетов £ = I0-2 (порядка 1 %) и ниже.