Численное моделирование пиролиза этана (2,1)-методом ре-шения жестких неавтономных задач Текст научной статьи по специальности «Математика»

УДК 519.622 Л.В. Кнауб, A.E. Новиков, Ю.А. Шитов

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПИРОЛИЗА ЭТАНА (2,1)-МЕТОДОМ РЕШЕНИЯ ЖЕСТКИХ

НЕАВТОНОМНЫХ ЗАДАЧ*

Построен L-устойчивый (2,1)-метод второго порядка точности для решения жестких неавтономных задач. На каждом шаге один раз вычисляется правая часть системы дифференциальных уравнений и осуществляется декомпозиция матрицы Якоби, два раза выполняется обратный ход в методе Гаусса. Эффективность метода интегрирования продемонстрирована на примере расчета пиролиза этана.

Ключевые слова: (т,к)-методы, оценка ошибки, контроль, точность вычисления, пиролиз этана.

L.V. Knaub, A.Ye. Novikov, Yu.A. Shitov

ETHANE (2,1) PYROLYSIS NUMERICAL MODELING BY THE TECHNIQUE OF STRICT NON-AUTONOMOUS

TASKS SOLVING

L-stable (2,1)-second precision order method for strict non-autonomous tasks solving is made. At every step differential equations system right part is calculated once and Jacobi matrix decomposition is also made; the Gauss technique return motion is completed twice. The efficiency of the integration technique is demonstrated on the ethane pyrolysis calculation example.

Key words: (m,k)-techniques, error estimation, control, calculations exactness, ethane pyrolysis.

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

В некоторых случаях расчеты требуется проводить с так называемой инженерной точностью - порядка 1%. Это связано с тем, что измерение констант, входящих в правую часть системы дифференциальных уравнений, часто проводится достаточно грубо. Иногда такая точность расчетов является удовлетворительной с точки зрения поставленной цели. Порядок аппроксимации численной схемы следует сочетать с требуемой точностью расчетов [1].

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

Здесь построен L-устойчивый (2,1)-метод второго порядка точности, в котором допускается замораживание как численной, так и аналитической матрицы Якоби. Приведены результаты численного моделирования пиролиза этана, подтверждающие эффективность алгоритма интегрирования.

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

У' = /(иУ\ У{10) = Уо> (1)

где у и 1 - вещественные М-мерные вектор-функции, I - независимая переменная.

Рассмотрим (2,1)-метод вида [4, 5]

'Работа поддержана грантами РФФИ №08-01-00621 и Президента НШ-3431.2008.9.

Уп+1=Уп+РА+Р2к2>

АЛ1 = КА*п + РКу„), АЛ = к1

(2)

с коэффициентами

/?, =б/ = 1- 0.5 л/2, р2 =0.5^2, = 0.5.

(3)

Коэффициенты (3) приведены для наглядности, а ниже их выбор будет обоснован. Здесь к и к2 -стадии метода, Оп=Е-аЬпАп, Е - единичная матрица, Лп - шаг интегрирования, Ап - некоторая матрица, представимая в виде Ап=Гуп+ ЬпВп +0(Л2п), ?у1п=д\(Пуп)1ду - матрица Якоби системы (1), Вп - не зависящая от шага интегрирования произвольная матрица, а, р1 и р2 - числовые коэффициенты. Использование при построении метода матрицы Ап позволяет применять (2) с замораживанием как аналитической, так и численной матрицы Якоби [6]. Если матрица Якоби ?у,п-к вычислена к шагов назад, имеем

где 0<к<1ь, 1ь - максимальное число шагов с замороженной матрицей, Вп=Г'у,п+ Т^п- Если матрица Якоби вычисляется численно с шагом г=с)ЬП, где с - постоянные, то получим Ап=^у,п+ЛпВп+0(Л2п), где элементы Ьщ матрицы Вп определяются по формулам

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

Теперь перейдем к исследованию схемы (2). Ниже для сокращения записи будем полагать Л=Лп, если это не будет приводить к недоразумениям. Разлагая к1 и к2 в ряды Тейлора в окрестности точки уп до членов с Л3 включительно

где элементарные дифференциалы вычислены на приближенном решении уп. Разложение точного решения у(^+1) в ряд Тейлора в окрестности точки и до членов с Л3 включительно имеет вид

Ап=дЦ1п-к,уп-к)1ду=Гу,п-кЬпВп+0(Ь2п),

(4)

Ьщ1=0,5с^п,уп)1ду21, 1 < , < N.

и подставляя в первую формулу (2), получим

уп+1 =у„ + (а+а)¥и + Да + Рг)Л* + «(а+2 а +

+/г3[0.5 Д(а + а)Л"и +яДа +2а)/'Х, +

+«2(а+3 а +я(а+2а)^и/и]+^/?4).

(5)

хы=хо+¥+о-5^2[^+/;л+

+(/гз/б)[у;;+/д'+2/;;/+/;2/+/;/2]+о(/г4),

где элементарные дифференциалы вычислены на точном решении у(Ц. Запишем условия второго порядка точности схемы (2). При условии уп= у(Ъ), сравнивая (5) и (6) до членов с Л2 включительно, получим

р1+р2= 1, а(р1+2р2) = 0.5, /? = 0.5. (7)

Теперь исследуем устойчивость схемы (2). Для этого применим ее к скалярному тестовому уравнению у'=Ку, у(0)=уо, [1е(Л)<0. Получим уп+1=0(х)уп, где х=ЛЛ, а функция устойчивости 0(х) с применением (7) имеет вид 0(х)=[1+(1-2а)х+а(а-р1)х2]/(1-ах)2. Отсюда условие ^-устойчивости (2) записывается в виде р1=а. Учитывая (7), из соотношения р1=а получим а2-2а+0,5=0. С использованием разложений стадий в ряды Тейлора локальную ошибку 5п метода (2) можно записать в виде

8п=(а- 1/3)й3/;7 + А3[^Л'+\/;,г +

Тогда согласно [7] для контроля точности (2) можно использовать оценку ошибки вида £п=(а-1/3)Ь%/+0(Ь3). Так как к2-к1= аЛ2/,у,п/п+0(Л3), то £п с точностью до членов 0(Л3) можно оценить по приближенной формуле

£(/'п)=а-1(а-1/3)Оп1-п[к2-к1], <<2.

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

1Нл)11^^/|я-1/3|, (9)

где е- требуемая точность интегрирования, \\\\- некоторая норма в ^, а целочисленная переменная /п выбирается наименьшей, при которой выполняется неравенство (9). Отметим одну важную особенность построенной оценки ошибки (9). Схема (2) ^-устойчивая, то есть для ее функции устойчивости 0(х) имеет место соотношение 0(х)^0 при х ^ -«. Так как для точного решения у(^+1)=ехр(х)у(п задачи у'=Ку, Яе(Л)<0 выполняется аналогичное свойство, то естественным будет требование стремления к нулю оценки ошибки при х ^ -«. Однако для е(1) это свойство не выполняется - данная оценка ведет себя А-устойчивым образом. С целью исправления асимптотического поведения ошибки вместо е(1) введена оценка е(/п), 1^ /п ^ 2. В этом случае поведение оценки ошибки при /п=2 будет согласовано с поведением точного решения тестовой задачи при х ^ -«. Подчеркнем, что в смысле главного члена оценки е(1) и е(2) совпадают. Использование е(]п) фактически не приводит к увеличению вычислительных затрат. Это связано с тем, что е/) при /п=2 проверяется только в том случае, если оно нарушено при /„=1. Такая ситуация встречается достаточно редко, в основном при быстром росте величины шага интегрирования.

Уравнение а2-2а+0,5=0 имеет два корня а1=1-0^2 и а2=1+0,5^2. Выберем а=а1, потому что в этом случае меньше коэффициент в главном члене (а-Ш^/у локальной ошибки (8). Теперь коэффициенты схемы (2) определяются однозначно и совпадают с (3), а соотношение а/|а-1/3| в неравенстве (9) приблизительно равно 7. При решении конкретных задач с целью повышения надежности расчетов коэффициент перед е в неравенстве (9) полагался равным единице.

Прежде чем формулировать алгоритм интегрирования, приведем некоторые соображения по выбору величины шага численного дифференцирования и алгоритму замораживания матрицы Якоби. При численном вычислении матрицы Якоби шаг численного дифференцирования / 1<]^ выбирается по формуле Г=тах(10-14, 10-7|у/|). Решение жестких задач осуществляется с двойной точностью в силу плохой обусловленности матрицы Якоби системы (1). Постоянная величина 10-7 введена в формулу определения шага численного дифференцирования с целью выдвижения его на середину разрядной сетки. Теперь /-й столбец а)п матрицы Ап в формуле (2) вычисляется по следующей формуле:

ап=--------------------------------------------------------, 1<у<А/\

О

то есть для вычисления матрицы Ап требуется N вычислений правой части задачи (1). Попытка использования прежней матрицы йп осуществляется после каждого успешного шага интегрирования. Для того чтобы не испортить свойства устойчивости численной схемы, при замораживании матрицы йп величина шага интегрирования тоже остается постоянной. Размораживание матрицы происходит в следующих случаях: 1) нарушена точность расчетов, 2) число шагов с замороженной матрицей достигло заданного максимального числа Ь, 3) прогнозируемый шаг больше последнего успешного в Оь раз, 4) е(1)>е(2). Поясним последнюю причину размораживания более подробно. Из результатов расчетов следует, что е(2) вычисляется либо при резком увеличении шага интегрирования, либо при невыполнении точности расчетов при быстром изменении решения. Если матрица Якоби вычисляется на каждом шаге, то е(2) более точно отражает асимптотическое поведение ошибки, и ее использование позволяет избежать неоправданных возвратов при резком возрастании величины шага интегрирования. Однако, если применяется матрица Якоби, вычисленная несколько шагов назад, то оценка ошибки е(2)=||йп'1(к2- к1)|| может оказаться хуже е(1) за счет применения "испорченной" матрицы. Поэтому, если е(1)>е(2), а неравенство (9) выполняется при /п=2, то полученное приближение к решению принимается, но следующий шаг выполняется с перевычисленной матрицей Якоби.

Теперь сформулируем алгоритм интегрирования. Пусть признак I определяет способ вычисления матрицы Якоби - численная или аналитическая. Предположим, что приближение уп к решению у(п задачи (1) вычислено в точке и с шагом Ьп, при этом к шагов осуществлялось с замороженной матрицей Якоби. Учитывая, что е(/п)=0(Л2), 1</п<2, алгоритм интегрирования формулируется следующим образом.

Шаг 1. Если к=0, то в зависимости от признака I вычисляется матрица Якоби.

Шаг 2. Формируется матрица йп и выполняется ее декомпозиция.

Шаг 3. Вычисляются к1 и кг.

Шаг 4. Вычисляется оценка ошибки е(1) по формуле е(1)=||к2- к1||.

Шаг 5. Вычисляется величина ф по формуле ф21е(1)=е, где е - требуемая точность вычислений, ||-|| -некоторая норма в ^.

Шаг 6. Если ф<1, то есть требуемая точность не выполнена, то вычисляется е(2) по формуле е(2)=||й 1п(к2-к1)||. В противном случае значение е(2) полагается равным е(1).

Шаг 7. Вычисляется значение ф2 по формуле ф22е(2)=е.

Шаг 8. Если ф2<1, то Лп полагается равным ф2Лп. Если к=0, то возврат на шаг 2; в противном случае к полагается равным нулю и происходит возврат на шаг 1.

Шаг 9. Вычисляется приближение к решению в точке ^+1 по формуле (2) с коэффициентами (3).

Шаг 10. Вычисляется значение Лпеш по формуле ^леи=т\п(ф, ф2)Лп.

Шаг 11. Значение к полагается равным к+1.

Шаг 12. Если к<1ь, Ьпеш<ОьЬп и е(1)<е(2), то Лп+1 полагается равным Ьп. В противном случае к полагается равным нулю и Лп+1=Лпе». Здесь 1ь - максимальное число шагов с замороженной матрицей Якоби, Оь - коэффициент роста величины шага интегрирования.

Шаг 13. Выполняется следующий шаг интегрирования.

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

Расчеты проводились на 1ВМ РС А№1оп(Ьт)ХР 2000 с двойной точностью при 4=20 и Оь=2. Норма в неравенстве (9) задавалась формулой ||k2-k1||=max1<<N|ki2-ki1|/(|yin|+f)|, где I - номер компоненты, г - положительный параметр. Если по 1-й компоненте решения выполняется неравенство |у'п|<г, то контролируется абсолютная ошибка ге, в противном случае - относительная ошибка е. В расчетах параметр г выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой. Вычисления проводились с точностью е=10-2. Это связано с тем, что в построенном алгоритме применяется схема низкого порядка точности, и поэтому данным методом осуществлять расчеты с более высокой точностью нецелесообразно. Сравнение эффективности проводилось с известным методом Гира сИвос1е в реализации А. Хиндмарша [8]. Ниже через // и / обозначены соответственно суммарное число вычислений правой части

и количество обращений (декомпозиций) матрицы Якоби задачи (1), которые позволяют объективно оценить эффективность алгоритма интегрирования.

Пиролиз этана в отсутствии кислорода описывается небольшой последовательностью стадий. Механизм пиролиза этана неоднократно обсуждался в литературе. Здесь принята схема реакции, предложенная и исследованная в работе [9].

С2Я6->ся3+ся3, ся3 + С2Н6 ся4 + С2Я5,

С2Н5 -> С2Я4+Я, Я + С2Я6^Я2+С2Я5, С2Я5+С2Я5^С4Я10.

Константы скоростей стадий имеют вид: ^=1,34^10-5, fc=3,73^102, fe=3,69^103, А'4=3,66-105, fc=1,62^107. Обозначим концентрации реагентов следующим образом: ci=[C2He], С2=[СНз], сз=[СН4], С4=[С2Н], С5=[С2Н4], С6=[Н], С7=[Н], C8=[C4Hio]. Тогда соответствующая система состоит из 8 обыкновенных дифференциальных уравнений вида

с[ = —kxcx — k2c]c2 - k4cxc6, с2 = 2kxcx — k2cxc2, c\ = k2cxc2,

c4 = k2cxc2 - k3c4 + k4cxc6 - 2k5c4, (10)

^3^41 Cg ^3^4 k4cxc61

c1 = k4cxc6, c8 = k5c4.

Начальная концентрация этана С1=[С2Нб] равна 0,14, для остальных реагентов начальные концентрации равны нулю.

Расчеты осуществлялись с точностью £=10-2. Матрица Якоби задачи (10) вычислялась численно. Эффективность алгоритмов интегрирования оценивалась по числу вычислений правой части if и матрицы Якоби ij задачи (11) на интервале интегрирования. Численное решение осуществлялось на промежутке [0; 0,26] с начальным шагом h=10-5. Данная задача удовлетворяет "классическому" определению жесткости. В начале интервала интегрирования наблюдается переходный участок (сотые доли секунды), а затем происходит медленное установление.

Решение задачи (10) построенным алгоритмом удалось вычислить с затратами if=9 и j=4. Для метода Гира dlsode затраты следующие - if=37 и j=3. При более высокой точности расчетов программа dlsode эффективнее алгоритма интегрирования переменного шага (2), (3). Это является следствием низкого порядка численных формул (2), (3).

Таким образом, построенный алгоритм предназначен для расчетов с небольшой точностью - порядка 1% и ниже. В этом случае достигается его максимальная эффективность.

Литература

1. Об оптимизации параметров методов численного решения задачи Коши для обыкновенных дифференциальных уравнений / С.С. Артемьев [и др.] // Численные методы механики сплошной среды. -1984. - №2. - С. 5-14

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

3. Byrne G.D., Hindmarsh А.С. ODE solvers: a review of current and coming attractions // J. of Comput. Physics, 1987. - №70. - Р. 1-62.

4. Новиков Е.А., Шитов Ю.А., Шокин Ю.И. Одношаговые методы решения жестких систем // ДАН СССР. - 1988. - Т. 301. - №6. - С. 1310-1314.

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

6. Новиков Е.А., Шитов Ю.А. Алгоритм интегрирования жестких систем на основе (т,к)-метода второго порядка точности с численным вычислением матрицы Якоби: Препринт №20. - Красноярск: ВЦ СО АН СССР, 1988. - 23c.

7. Демидов Г.В., Новиков Е.А Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений // Числ. методы механики сплошной среды. - 1985. - Т. 16. - №1. - С. 27-42.

8. http://www.netlib.org/odepack/index.html.

9. Kukh D.M., Taylor J.E. Mathematical simulation of the oxygen ethane reaction// J. Chem. Kinet. - 1975. -V.8. - P.89-97.