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

раздел МАТЕМАТИКА

УДК 541.127 : 519.62 : 517.928

ПРИБЛИЖЕННЫЕ ФОРМУЛЫ РЕШЕНИЯ ПРЯМОЙ КИНЕТИЧЕСКОМ ЗАДАЧИ ДЛЯ ЖИДКОФАЗНОГО ИНГИБИРОВАННОГО ОКИСЛЕНИЯ УГЛЕВОДОРОДОВ

© А. В. Тропин *

Башкирский государственный университет Россия, Республика Башкортостан, г. Уфа, 450074, ул. Фрунзе, 32.

Тел./факс: + 7 (34 7) 273 6 7 78.

* E-mail: TropinA V@rambler.ru

Эта статья является второй в цикле трёх статей автора по данной тематике. Использование алгоритма, описанного в предыдущей работе «Алгоритм построения приближённых формул для решения прямой кинетической задачи на основе численного решения» иллюстрируется на известном примере. Подробно обсуждаются все этапы алгоритма: численное интегрирование, упрощение исходной задачи, получение локальных формул, сращивание их в глобальные, проверка малой зависимости погрешности последних от констант скоростей.

Ключевые слова: прямая кинетическая задача, жесткость, численное решение, точное решение, приближенные формулы, правило упрощения.

Рассматривается упрощённый (для краткости) вариант механизма жидкофазного цепного окисления органического субстрата КН молекулярным кислородом 02 в присутствии ингибитора 1пН с участием радикалов Я02, 1п [1].

RO 2 , k 0 » 10 8 моль /(л ■с)

RO 2 + RO 2 ^ , k 6 » 10 6 и л /(моль с)

InH + RO 2 In , k 7 » 10 4 л /( моль с)

RO 2 + In , k 8 » 10 9 л /( моль с)

In + In , k 9 » 10 8 л /( моль с) .

Здесь инициирование процесса представлено в виде формальной стадии с константой k0. Опущены продукты некоторых стадий. В соответствующих размерностях заданы ориентировочные значения констант скоростей.

Пусть начальные концентрации имеют следующие порядки [inH]0 = 10- 7 ■ моль / л,

[RO2 ]0 = 10 5 ■ моль / л, [In ]0 = 0 ■ моль / л

Введём обозначения для текущих и начальных концентраций

[InH] = у, [RO2 ] = х, [In ] = z, [InH]0 = y0, [RO2 ]0 = x0

В этих обозначениях, согласно [2], изменение концентраций во времени описывается задачей Коши (системой обыкновенных дифференциальных уравнений (СОДУ) химической кинетики с начальными условиями)

d = -k7 ух, y(0) = V

dx = k_ -2k,х2 -k-ух-k„xz, х(0) = хп dt 0 6 r 8 0 .

dz = knyx - k„xz - 2knz2, z(0) = 0 dt r 8 9

(1)

Точное решение задачи (1) выписать затруднительно. В то же время можно найти для него приближенные формулы, достаточно точные на [0,да). Для этого исходная задача Коши (1) предварительно интегрируется численно, например, яв-

ным методом Рунге-Кутта с переменным шагом (пропорциональным логарифму от времени), но с неявной схемой Эйлера для “быстрых” переменных. Теория жёстких СОДУ [3] и опыт автора показывает, что такое сочетание наиболее эффективно при численном решении задач Коши в химической кинетике.

Задаётся порог малости е=0.01. На основе численного решения выясняется на каких интервалах времени и какие члены в СОДУ задачи (1) могут быть отброшены так, чтобы возмущение решения, вносимое этим упрощением, было мало (=е). В работах [4,5] выяснено, что СОДУ в общем виде (За. г

1н = . ^ Д ]"]<а<,))

/ = 1

рационально упрощать по следующему правилу.

Если для некоторого /', при / из интервала времени [а, Ь) верно неравенство

r

2 У: Г W (a(t)) < e ■ max )t a(( £

J = 1 г,J J i, J J

то 1-е дифференциальное уравнение на [а, Ь) можно заменить на алгебраическое г

0 = 2 У- ■ (а(0)

1 1, ] ] .

] = 1

При этом, если для некоторого ] выполняется неравенство

g . ■w ,(a(t))

i, J J

g . ■ w . (a(t))

i, J J

< e ■ max g J

то слагаемое yij-wj(a(t)) в этом алгебраическом уравнении можно удалить.

Если для некоторых i и J при t е [a, b) выполняются оба неравенства

r

M g w aa > e ■ max £

j=1, j j J ‘, J J

g . ■ w . (a(t)) < e ■ max<

\l, J J I i

M g ,- W. (a(t))

J = 1 ', J J

то на [a, b) в i-м дифференциальном уравнении

/ max(a. (t)) > ■ max(a. (t))

t 1 I t 1

r

6

раздел МАТЕМАТИКА

слагаемое а(/)) можно удалить. Когда аі

практически не изменяется, получается уравнение вида dai/dt=0.

Применяя вышеизложенное правило к значениям численного решения в узловых точках, выявлено 3 интервала времени [0, 0.1), [0.1, 10), [10,да), на которых исходная задача Коши (1) упрощается до достаточно простых задач Коши, допускающих точное аналитическое интегрирование.

Так, на интервале [0, 0.1) исходная задача упрощается до следующей

| - °- у(0) - >0. 1 - 0 х(0 = -0І

dz = knyx - kaxz, z(0) = 0 dt r 8 .

Её решение в системе Maple легко находится:

У = Уг

z = (1 - exp^- k8x0t j j ■ y0k7 / k8

*0’ ~ ~o’ ^ '^0jJ 7rr 8.

На интервале [0.1, 10) исходная задача

упрощается до системы

d=0, у(0)=У0; d=ko-2k6x2-k7yx, x(0)=xo;

0 = k7 yx - k^xz.

Находя г из алгебраического уравнения, подставляя его в дифференциальное, программа находит решение системы:

1 \ - к7 + 1апЬ

2k

t■ \l2k0k6 + k72 У02 +

+ 1ln

2

2k0k6 + k72 y02 + k7 y0 + 2 x0k6

2k0k6 + k72У02 - k7y

2x k

0 6

2k0k6 + k72 yo2

y = y0 z = y^ki/ k

0 7 8

На интервале [10,да) упрощённая задача имеет

вид

^ - -к7ух, >(0) - >0 0 - ^ - 2к6х2 - ух - к^хх .

0 - к,ух - к0хг.

/ 8

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

к7(у-у0)+1|(к7у)2 + 2к0к6 ^(к7у0)2+2к0к6 ^2к0к6 • 1п(У/у0)-

^ ^+^+(к7у^2к0к^1 1+((^>п )/(2^^)

1/2k0k6 ’ AT ’ «~7У)'(2~0"6^ Д1 ' АТ ’ '.("7У0)'("''0"'6

■J(k7 у)2+2k0k6 - k7yj /(2k6); z = y V V

-k0k7t

■-7^ -0 6 ^ у 6’

Интересно заметить, что при соотношении параметров (к7у0)2 >> 2к6к0 выражение для у можно приближённо упростить до явного вида.

Полученные выше формулы решения трёх упрощённых задач Коши легко срастить (скомби-

нировать), то есть выписать глобальные формулы, достаточно точные на всём [0,да) для исходной задачи. Причём опыт показывает, что при сращивании относительные погрешности глобальных формул всегда оказываются намного меньше, чем погрешности исходных локальных формул.

Сопоставим формулы для у на всех трёх временных интервалах. Рассмотрим неявное уравнение для у на интервале [10,да). При t<10 модуль правой части k0k7t < 10-3. Значит, и левая часть уравнения тоже должна быть мала по модулю. В Maple легко проверить, что это возможно лишь в случае у~у0. Следовательно, неявное уравнение для у на интервале [10,да) будет достаточно точно и на интервалах [0, 0.1), [0.1, 10).

Сопоставим формулы для х на всех трёх временных интервалах. Рассмотрим выражение для х на [0.1, 10). При t>10 аргумент гиперболического тангенса велик, значит, tanh(...) близок 1, а тогда это выражение почти совпадает с формулой для х на интервале [10,да). В системе Maple также достаточно легко проверить, что при t<0.1 то же выражение для х близко к х0. Для этого достаточно воспользоваться эквивалентностью бесконечно малых из теории пределов.

Совсем легко проверить, что выражение для z на [0, 0.1) близко к k7y/k8 при t>0.1, так как отрицательная экспонента мала.

Таким образом, глобальные формулы могут выглядеть, например, так:

k7(у-у0 W(k7у)2 + 2k0k6-\(k7Уо)2 + 2k0k6 +V2kok6 ■ ln(y7Уо)-

F]) = -k0k7t ’

7^' 0 6 V 7-^' 0 6

Ы2k0k6 lnH141 + ((k7y)/(2k0k6))2 'I/f1 + М(^У0)/(2^6).

■tanhl t

x=^{- k7У 4

\2.kok^ + kT2 y2

h( *■ J

2k0k6 + k72 y2 +

l2kr,k*+ki2 у 2 + kiy+2xA

+k72 y 2 -

k7y - 2x0k6

г = ^1 - ехр|^- к^х^^ ■ у ■ к7 / к8 .

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

Для выяснения области значений коэффициентов к0, к6, к7, к8, к9, при которых формулы достаточно точны, были проведены численные эксперименты.

В начале коэффициенты варьировались по отдельности множителями в интервале [0.1, 10). Вычислялись предельные относительные погрешности формул при изменённых значениях констант. В качестве точного решения использовалось численное. Результаты эксперимента представлены на рис.1 в виде графиков зависимости корня из суммы квадратов предельных относительных погрешностей нормированных концентраций от варьирующего множителя для каждой из констант.

/

Рис.1. Зависимости предельной относительной погрешности формул от варьирующего множителя при варьировании констант скоростей.

Из рисунка 1 видно, что предельная относительная погрешность приближённых формул быстро нарастает при увеличении к7 и уменьшении к8. Важно заметить, что само решение существенно изменяется в процессе варьирования.

Для выяснения причин роста погрешности на рис. 2 совмещены графики приближённых аналитических формул и точного решения при увеличении к7 в 10 раз. Ось абсцисс логарифмическая, соответствует времени. По оси ординат откладываются концентрации, нормированные на максимальные.

Рис.2. Сравнение численного решения исходной задачи (сплошная линия) и приближённого по формулам (пунктир).

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

(приближённое) сильно расходится с г(точным). Значит, эту компоненту в приближённых аналитических формулах желательно подправить. Это удаётся сделать по той же методике, оставляя слагаемое -2к9г2 третьем уравнении упрощённой подсистемы на интервале [0, 0.1). При этом выражение для г несколько усложнится. Менее чувствительным к увеличению к7 и уменьшению к8 будут глобальные формулы с подправленным выражением для х и существенно изменённым для г.

к7 (У -у0 М2к0к6 + (к7у)2 -т/2к0к6 + (к7У0)2 +Т2к0к6^ 1п(у/у0 )-

^2к0к6 ■ 1п[(1^ 1+((к7у)/(2к0к6))2 )^+^ 1+((к7у0)/(2к0к6))2)) = -к0 к7;

+ - 1п

2

2к0к6 + к72 у02) + к7 у0 + 2х0к6

2к0к6 + к72 V)- к7 у0 - 2х0к6

г = |- к8 х + д/ 8к7 к9 хУ + к« 2x2 ■‘ап4 ^л/8к7кохпУп + к« 2хп 2 +

8

2\ Г9 0^0 8 0

+ 121п {Гк7к9х0у0 + к82х02 + к8х0 )/( V 8к7к9х0у0 + к82х02 - к8х0

Для этих формул относительные погрешности значительно меньше. Это ясно из рис.3 - аналога рис.1 для подправленных формул.

0.025

0.02 /4-7

0.015 \А'<5 у кб

0.01

0.005 к9

4-7 к8

х=2Т Г к7 у Ч 2кок6+к72У2 ■ 1ап^ *■ '12к"к* + ^ 2^ 2 +

0 6 7 ^0

0.1 0.2 0.5 1 2 5 10

Рис.3. Зависимость предельной относительной погрешности уточнённых формул от варьирующего множителя при варьировании констант скоростей.

Во втором эксперименте все коэффициенты к0, к6, к7, к8, к9 варьировались одновременно. Исходные значения независимо умножались либо на 0.1, либо на 10. Перебор всех вариантов умножения показал: из 32 вариантов варьирования 21 раз получилась погрешность меньше 0.1.

В третьем опыте параметры варьировались так же, но варьирующие множители для каждого из параметров выбирались случайно из отрезка [0.1, 10]. В 120 случаев из 200 относительная погреш-ность оставалась меньше 0.005.

Из этих экспериментов следует, что полученные формулы достаточно точно описывают кинетику процесса не только для фиксированного набора значений констант к0, к6, к7, к8, к9 , но и для значений близких порядков. В отличии от численного решения, эти формулы не являются просто количественной аппроксимацией информации о процессе. Поэтому, на наш взгляд, более естественно использовать именно их при моделировании химических процессов. Например, используя их при решении обратной кинетической задачи, можно резко сократить объем оптимизационных вычислений.

ЛИТЕРАТУРА

1. Эмануэль Н. М., Денисов Е. Т., Майзус З. К. Цепные реакции окисления углеводородов в жидкой фазе. М.: Наука, 1965. -375 с.

Вант Гофф Я. Г. Очерки по химической динамике. Л.: ОНТИ ХИМТЕОРЕТ, 1936. -178 с.

К. Деккер, Я. Вервер. Устойчивость методов Рунге-Кутты для жёстких нелинейных дифференциальных уравнений. М.: Мир, 1988. -336 с.

Тропин А. В., Спивак С. И. // Кинетика и катализ. 1995. Т.36. №5, С. 658-664.

Тропин А. В. Редукция математических моделей механизмов цепных реакций: дис. ... канд. физ.-мат. наук. Уфа, 1998. -104 с.

Поступила в редакцию 11.04.2007 г.