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

УДК 541.127:519.62:517.928

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

© А. В. Тропин

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

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

E-mail: TropinAV@rambler.ru

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

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

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

Рассматривается сложный цепной химический процесс [1]

■ O3 OH■ + R ■ , k1 » 10-2 c_1;

O3 , k2 » 0.0162 c_1;

OH■ R ■ , k3 » 9 ■ 105 c-1;

r ' ro2 , k4 » 2.4 407 c_1;

O3 + r ro■ , k5 » 7 408 л /(моль ■c);

O3 + RO2 ■ RO■ , k6 » 1.2404 л /(моль ■ c);

RO■ R ■ , k7 » 9 105 c_1;

R ■ + R ■ , k8 » 2.1109 л/(моль ■c);

RO2 ■ + R ■ , k9 » 1.5 109 л /(моль ■c);

RO2 ■ + RO2 ■ , k0 » 2 406 л /(моль ■c).

Пусть значения начальных концентраций имеют следующие порядки (моль/л):

- 4

[O3]0

: 0, [OH ]0 = о, [R ]0 = 0.

10 4, \яо2]0 = 0, [ЛО']0

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

[Оз] = и, [ЯОг] = [КО] = [ОЯ ] = х, [К ] = V, Оз]0 = и0.

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

Лы/ Л = -к^и - А^и - к и - к^иу, и(0) = ^ ;

2

dу|dt = k4V-к^иу-к^уV-2^у , у(0) = 0; сЬ/ dt = к и + к^иу - к/ z, z(0) = 0 ■ dх|dt = к^и - кдХ, х(0) = 0 ■

dv|dt = к.и + к^х-к^-ки + кп2-2к^ -к,у, v(0) = 0 13 4 3/оУ

Очевидно, что точное решение этой задачи Коши выписать невозможно. Однако можно найти для него приближенные аналитические формулы, достаточно точные на всём интервале [0,да). Для этого исходная задача Коши сначала интегрируется численно каким-либо методом, нечувствительным к типичной жёсткости СОДУ химической кинетики. Подробности этого рассмотрены в предыдущих двух статьях, а проблема жёсткости рассмат-

ривается полно в [3].

Задаётся порог малости є = 0.01. На основе численного решения программа выясняет, на каких интервалах времени и какие члены в СОДУ задачи Коши могут быть отброшены так, чтобы возмущение решения, вносимое этим упрощением было мало (порядка є). В работах [4, 5] выяснено оптимальное правило упрощения СОДУ. В предшествующих данного цикла [6, 7] оно дважды строго сформулировано. По этому правилу принимает решение программа. При применении вышеупомянутого правила к значениям численного решения в узловых точках выявляется три интервала времени [0, 10-4), [10-4, 1), [1,®), на которых исходная задача Коши может быть упрощена до задач, допускающих точное аналитическое интегрирование.

Итак, на интервале времени [0, 10-4) исходную задачу можно упростить до

— = 0, u(0) = u„ dt 0

— = k.uv - k_z, dt 57

dy ,

— = k4v,

dt z(0) = 0

y(0) = 0

— = ku-k,x, x(0) = 0 ■ — = k.u + k„x-k.v, v(0) = 0. dt 1 3 ■ dt 1 3 4

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

имеет следующий явный вид:

u(t) = u

0

y(t) = k1uQ

2 + k4 exp(-k3t) + (k4 - 2k3) exp(-k4t) (k4 + 2k3)

k3(k4 -k3)

k4(k4 -k3)

k3k4

z(t) = kku 2| +

(k4 - 2k3)exp(-k4t) (2k3 - k7)exp(-k7t)

I"5"0 [k/7 ' k - k,)(k7 - k,) + k4(k4 - k3)(k4 - k7) “ k7(k4 - k7)(k3 - k7)

exp(-k3t)

4 - кз)(к/ - кз

х(t) = к1и01 - exp(-kзt))/ кз;

( 2 exp(-k3t) (к4 - 2k3)exp(-k4t)

v(t) = к,иЛ---------------------------—

и 10[к4 (к4 -к3) к4(к4 -к3)

На интервале времени [10-4, 1) СОДУ исходной задачи Коши сокращается до подсистемы дифференциально-алгебраических уравнений, поэтому

— = 0, u(0) = u„ dt 0

— = k.v - k,.uy - 2k„y2, y(0) = 0 dt 4 6 у 0* ’

0 = k.uv + ku-k„z ■ 0 = k.u -k~x ■

5 6 1’ 1 3 ’

0 = k^ + k3x - k^v - ku + ki z

Здесь три последних уравнения - алгебраические. Программа находит корни алгебраической подсистемы, подставляет эти корни в дифференциальное уравнение и выдаёт

— = 0, u(0) = u„ ■ dy = 2ku -2k„y2, y(0) = 0 ■ dt 0 ■ dt 1 0 ■

2 2

z = (Ik.k.u + k.k,u y + k.k^uy)/(k.k„) ■ x = k.u/k0 ■

15 5 6 4 6 4 7 ’ 1 3 ’

v = u ■ (2k1 + k6 y) / k4.

Далее Maple легко находит решение сокращённой задачи Коши. Окончательно на интервале времени [10-4, 1) получается следующее решение:

u(t) = u0 ■ y(t) = (k1u0/k0)1/2tanh|2t^(k0k1M0)1/2' ■ z(t) = u0 ■ (2k1k5u0+k5k6u0y(t)+k4k6y(t))/(k4k7) ■ x(t) = k1u0 / ■ v(t) = u0 ■ (2k1 + k6y(t))/ k4

На интервале времени [1, да) СОДУ исходной задачи Коши можно сократить до подсистемы дифференциально-алгебраических уравнений в виде:

— = -k.u -k„.u -kruy, u(0) = u„ ■ dt 1 2 6 0

2

0 = k4v - k^uy - 2^y ■ 0 = k^uy - kiZ ■ 0 = k^u - k^x ■

0 = k^ + k3x - k4v + ki z.

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

u(t) = u0(k1 + k2)2 exp(-t(k1 + k2)^^k1 + k2 + ^(k^ /kQ)1/2 <1- exp(-t(^ + k2)/2) y(t) = (k^u(t)/^)1/2 ■ z(t) = k6«(t)/ki tk.ju.tt')/^)1/2 ■ x(t) = k^u(t)/k3 ■ v(t) = u(t)(2k.1 + k^k^t)/k^^)/k4 .

Полученные выше три решения хорошо сращиваются. Для этого рассмотрим поочерёдно выражения для каждой из концентраций на трёх интервалах.

Рассмотрим выражения для u(t). На интервалах времени [0, 10-4) и [10-4, 1) они неизменны u(t) = и0. Т акже легко заметить, что выражение для u(t) интервала [1, да) близко к u0 при t из первых двух интервалов. Поэтому естественно положить, что для всех te [0, да)

u(t) = M0(k1 + k2)2exp(-t(k1 + k2)^+ k2 + ^(k^/kQ)1/2■ (1-exp(-t(k1 + k2)/2))'

Срастим выражения дляy(t). Пусть te [0, 10-4). Тогда первое слагаемое в выражении для y(t) решения на интервале [0, 10-4) и решение на интервале [10-4, 1) близки, так как tanh(a) ~ а при а ^ 0. Значит первое слагаемое в выражении для y(t) 1-го интервала можно заменить на решение на втором интервале. И такое решение останется приблизительно верным на втором интервале, потому что остальные слагаемые малы при te[10-4, 1). Оно также будет верно на третьем интервале, т.к. tanh(a) ~ 1 при больших а. По понятным причинам некоторые u0 можно (и следует) заменить на u(t). Перемещая экспоненты в числитель

дробей и упрощая, можно получить выражение для у(/) пригодное на всём интервале [0, да):

кл ехр(-к^) (к4 - 2к3) ехр(-кд Г) (кл + 2к3)

y(t) = НТТГТТГ+

4‘) у '4 1 —3-

к3(к4- к3) к4(к4- к3) к3к4

+ 1апЪ^ 2г(к0к1^0)1/2 ^ ■ (к^ (0/ к0)1/2.

Оценивая порядки слагаемых в этом выражении при / из различных временных интервалов, легко проверить, что на каждом из интервалов оно близко к соответствующим локальным выражениям для у(г).

Рассмотрим выражения для 2({). Рассуждая аналогично, имеем для ?е [0,да)

z(t) = k.k,un 21 —2—I--

150 lkA (k

exp(-k3t)

(k4 - 2k,)exp(-^t) (2k, -kn)exp(-k7t)

3 7

k,)(k7 - k,) k4(k4 - k3)(k4 - k7) k7(k, - k7)(k, - k7)

7 4 7 3

7

Действительно, с помощью Maple можно проверить, что на временном интервале [0, 10-4) порядок второго слагаемого значительно меньше порядка первого слагаемого. Поэтому при te [0, 10-4) получаем в точности решение для z(t) на этом интервале. На интервале [10-4, 1) становятся малы слагаемые с экспонентами (внутри больших скобок) и получаем z(t) на этом интервале. На интервале [1, да) остаётся только первое слагаемое, которое мало отличается от решения для z(t) на этом интервале (так как u(t) << u0 и y(t) ~ (k1u(t)/k0)1/2).

Легко проверить, что x(t) сращивается следующим образом x(t) = A'1i/(t)(1 - exp^k^t))/ k^

так как при te[10-4, да) и k3~9-105, отрицательная экспонента ничтожно мала.

Аналогично получаем сращивание для v(t). Для te [0,да)

v(t) = u(t)

2k^ + k6 y(t) k^ exp(-k3t) + k^ (k4 - 2k3) exp(-k4t)

k4+k9 y(t)

k4(k4 - k3)

Оно тоже близко к локальным решениям на соответствующих интервалах.

Окончательно имеем следующие приближённые формулы для решения:

u(t) = -

^0 (k1 + k2) exp(-t(k1 + k2))

Л/2

k1 + k2 + k6(k1M0 / k0) (1 - exp(-t(k^ + k^)/2)):

y (t) = ^u (t)

k 4 exp( - k 31) (k 4 - 2 k 3 ) exp( - k 41) (k 4 + 2 k 3)

k 3 k 4

^ I k3(k4 - k3) k4(k4 - k3)

+ tanh 12t (k 0 k1 u 0)1/2 (k1 u (t)/ k 0)1/2,

2 exp(-k3t) (k. - 2k3)exp(-k4t) (2k3 - k7)exp(-k7t)']

z(t) = k,k, ^ 2I-------------------------------------------------------------------+-3-4-3-4-3-7-Ll +

150 Ik4k7 (k4 -k3)(k7 -k3) k4(k4 -k3)(k4 -4) k7(kA -k7)(k3 -k7) J

x(t) = k1u(t)(1 - exp(-k3?))/ k^

v(t) = u(^)

2k1 + k6 y(t) k1k4exp(-k3t) + k1(k4 - 2k3)exp(-k4t)

k4 + k9 y(t)

k4(k4 - k3)

+

+

k4k-u(t) y(t)

+

K(k4 - k5u(t))

k4k6u(t) y(t)

+

k7(k, -kjiit))

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

Поскольку все приближённые формулы для решения получились явными, то несложно визуализировать изменение их относительных погрешностей. Для этого в исходной СОДУ делается замена переменных

u=u(t)+Au(t), y=y(t)+Ay(t), z=z(t)+Az(t), x=x(t)+Ax(t), v=v(t)+Av(t).

Здесь u, y, z, x, v - компоненты неизвестного точного решения исходной задачи Коши, u(t), y(t), z(t), x(t), v(t) - полученные приближённые аналитические формулы, Au(t), Ay(t), Az(t), Ax(t), Av(t) -искомые абсолютные погрешности (от них легко перейти к относительным). После замены и упрощения получается СОДУ относительно неизвестных функций Au(t), Ay(t), Az(t), Ax(t), Av(t).

Подстановкой t = 0 в приближённые формулы легко проверяется, что в нуле они имеют нулевую погрешность. Поэтому начальные условия для абсолютных погрешностей будут нулевыми: Au(0) = 0, Ay(0) = 0, Az(0) = 0, Ax(0) = 0, Av(0) = 0.

Задача Коши строилась автоматически по программе в Maple. Численное решение полученной задачи, нормированное на максимальные концентрации, представлено на рис.1. Ось абсцисс логарифмическая, соответствует времени. По оси ординат откладываются относительные погрешности 5u, 5y, 5z, 5x, 5v для каждой из приближённых формул u(t), y(t), z(t), x(t), v(t) соответственно.

0.0025

0.002

0.0015

0.001

0.0005

0

■0.0005

iff

5у \5z

8^,5Л',5д;,8г б./.йХ \

1е-9 1е-8 1е-7 1е-6 1е-5 1е4 1е-3 1е-2 0.1 1 10 100

Рис. 1. График изменения во времени относительных погрешностей приближённых формул

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

Для выяснения области значений коэффициентов к1, к2, к3, к4, к5, кб, к7, к8, к9, к0, при которых приближённые аналитические формулы достаточно хорошо аппроксимируют точные решения исходной системы, были проведены численные эксперименты. Так же, как в примере, описанном в предыдущей статье, вначале коэффициенты варьировались по отдельности. Варьирование производилось множителями в интервале [0.1, 10). Предельные относительные погрешности приближённых формул оценивались исходя из решения вышеописанной задачи Коши с неизвестными функциями Дм(?), ДХ0, Д^(0, ЛХ0, Ду(0. При этом искали максимумы относительных погрешностей. Интересно отметить, что использование такого подхода к поиску погрешности (в отличие от примера предыдущей статьи, где каждый раз искалось численное решение исходной задачи Коши при новых параметрах) показывает, что СОДУ относительно абсолютных погрешностей значительно менее жёсткая, по сравнению с исходной СОДУ. Из этого следует меньшая трудоёмкость численного интегрирования.

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

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

Из рисунка видно, что предельная относительная погрешность полученных приближённых аналитических формул быстро нарастает при увеличении кб, к2 и уменьшении к0. Это может быть поводом для дальнейшего уточнения приближённых формул по тому же алгоритму. Для этого полезно сначала выяснить, на каких временах возникает большая погрешность, а потом попытаться исправить формулы, подобрав подсистему дифференциально-алгебраических уравнений с меньшим числом выброшенных членов. Возможно, где-то следует оставить члены с константами кб, к2 и к0. Пример такого подбора уже приводился во второй работе цикла, поэтому повторяться не будем.

Выводы

Исследование предложенного в данном цикле статей алгоритма на многих других примерах показывает его эффективность. Во всех рассмотренных других примерах всегда удавалось подбирать приближённые аналитические формулы для решения задачи Коши. Даже на первый взгляд в «безнадежных» случаях с дробными стехиометрическими коэффициентами (с наличием интервалов времени, на которых подсистемы не имеют точного аналитического решения) удавалось подобрать такие формулы. Естественно, что полученные формулы оставались достаточно точными при варьировании значений всех её параметров в достаточно широких диапазонах, несмотря на то, что само точное решение задачи Коши при этом существенно изменялось. А если при каком-то наборе значений параметров приближённые аналитические формулы становилось непригодным, то можно было их откорректировать или подобрать другие. Для этого достаточно было применить данный алгоритм заново, при новом наборе значений параметров.

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

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

задачи. Совершенно очевидно, что при этом подходе к последней резко сокращается объем оптимизационных вычислений, так как не требуется многократно решать прямую кинетическую задачу. Обратная кинетическая задача сводится к задаче нелинейной регрессии.

Разумно предположить, что данный алгоритм поиска приближённых формул для решений задач Коши применим не только в химической кинетике. Потенциальными областями применения алгоритма могут быть модели со сложными нелинейными СОДУ в теории массового обслуживания. Для таких СОДУ даже программа автора (реализующая алгоритм) абсолютно не требует переделок. Основное требование к СОДУ: их правые части должны быть суммой однотипных членов. Как показывает практика автора, это могут быть не только квадратичные члены, но и члены с дробными показателями, например. Другое естественное требование - устойчивость решения задачи Коши, иначе бессмысленно искать для него приближённые формулы.

Нельзя не отметить, что, как показывает практика, для задачи Коши с жёсткой нелинейной СОДУ есть во много раз больше шансов получить приближённые формулы для решения, чем для задачи с нежёсткой СОДУ. В значительной степени именно благодаря жесткости СОДУ удаётся удачно её упрощать.

ЛИТЕРАТУРА

1. Семёнов Н. Н. Цепные реакции. Л.: Госхимтехиздат, 1934.

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

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

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

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

6. Тропин А. В. // Вестник Башкирск. ун-та. 2007. Т.12. №3. С.15-17.

7. Тропин А. В. // Вестник Башкирск. ун-та. 2007. Т.12. №4. С.5-7.

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