автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Компьютерная реализация задач редукции систем дифференциальных уравнений химической кинетики

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

На правах рукописи

Вайман Александр Маркович

КОМПЬЮТЕРНАЯ РЕАЛИЗАЦИЯ ЗАДАЧ РЕДУКЦИИ СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ХИМИЧЕСКОЙ КИНЕТИКИ

05.13.16 - Применение вычислительной техники, математического моделирования и математических методов в научных исследованиях.

Диссертация на соискание учёной степени кандидата технических наук.

Научный руководитель : доктор физико-математических наук, профессор Спивак С.И.

Уфа - 1999

Введение 1

1. Литературный обзор. 4

1.1. Конструирование систем дифференциальных уравнений 4 химической кинетики (СДУ ХК)

1.2. Методы решения прямых задач для СДУ ХК. Основные трудности 7

1.3. Асимптотические методы в применении к СДУ ХК 12 1.3.1 Приближение безразмерной СДУ ХК цепочкой подсистем 13 дифференциально-алгебраических уравнений

1.3.1.1. Анализ точности приближения 13

1.4 Алгоритмы численного интегрирования жёстких СДУ ХК 20

1.4.1 Численные схемы 20

1.4.2 Компьютерные алгоритмы. 26

2. Постановка задачи компьютерной реализации. 27

2.1 Методы редукции СДУ ХК. 27

2.2 Техническое задание. 28

2.2.1 Назначение разработки 28

2.2.2 Требования к программному продукту 29

3. Математическое обеспечение задач редукции СДУ ХК 31

3.1 Описание математической модели 31

3.2 Алгоритм 33 3.3. Программная реализация 34

3.3.1 Общие сведения 34

3.3.2. Функциональное назначение. 34

3.3.3 Используемые технические средства. 35

3.3.4 Вызов и загрузка. 35

3.3.5 Входные данные 35

3.3.6 Выходные данные 37 3.4 Спецификации модулей программы. 39

4. Реализация программного продукта при построении кинетических моделей конкретных реакций. 47

4.1 Алкилирование ксиленола метанолом. 47

4.2 Парофазное окисление дурола в пиромеллитовый диангидрид 79

4.3 Реакция жидкофазного цепного автоокисления углеводородов 126

5. Выводы 128

6. Литература 129

7. Приложение 1 (текст программы) 135

8. Приложение 2 (тестирование) 153

9. Приложение 3 (справки об использовании результатов диссертационной работы) 169

Введение

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

Ранее был предложен алгоритм редукции систем дифференциальных уравнений химической кинетики, основанный на идее разделения времён. Была доказана его математическая корректность.

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

Цель работы.

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

Научная новизна работы.

Создано компьютерное обеспечение задач редукции систем дифференциальных уравнений химической кинетики большой размерности.

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

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

- анкетирование ксиленола метанолом.

- реакция жидкофазного цепного автоокисления углеводородов

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

Практическая ценность работы.

Высокий уровень сервиса программного продукта позволяет использовать его конкретному пользователю-химику при построении кинетических моделей широкого спектра кинетики сложных химических реакций : гетерогенного, гомогенного, ферментативного катализа, жидкофазного окисления и т.д. В настоящее время программа является активно используемым инструментом при математическом моделировании кинетики и термодинамики сложных реакций в Институте нефтехимии и катализа АН РБ и УНЦ РАН, Институте органической химии УНЦ РАН, на кафедре физической химии и химической экологии БГУ, кафедре математического моделирования и информационных технологий Бирского госпединститута.

Апробация работы.

Результаты работы обсуждены на конференциях : XIII Международная конференция по химическим реакторам Химреактор-13 (Новосибирск, 1996), международная конференция «Математические методы в химии и химической технологии» (Тула, 1996), всероссийская научная конференция «Физика конденсированного состояния» (Стерлитамак, 1997), третья международная конференция «Дифференциальные уравнения и их приложения», (Саранск, 1998 ), конференция «Математические методы в химии и технологиях (ММХ-11)» (Владимир, 1998), XIV международная

конференция по химическим реакторам (Томск, 1998.), the third International conference on Unsteady-State Processes in Catalysis, (Russia St. Petersburg 1998), международная научная конференция «Оптимизация численных методов», посвященной 90-летию рождения Сергея Львовича Соболева (1908-1998) ( Уфа, 1998), конференция «Современные проблемы естествознания на стыках наук» (Уфа, 1998).

Работа обсуждалась на научных семинарах кафедры математического моделирования и кафедры дифференциальных уравнений БГУ, отдела вычислительной математики Института математики УНЦ РАН, лаборатории математической химии ИНК АН РБ и УНЦ РАН, кафедры вычислительной математики и кибернетики УГАТУ, кафедры математического моделирования и информационных технологий БирГПИ.

Публикации.

По теме диссертации опубликовано 3 статьи [35 - 37] и 10 тезисов Международных и Всероссийских конференций [38 - 47].

1. Литературный обзор 1.1. Конструирование систем дифференциальных уравнений химической кинетики (СДУ ХК) Основным понятием химической кинетики является понятие скорости реакции. Для элементарной реакции это - число элементарных актов химического превращения, происходящих в единице реакционного объёма или на единице реакционной поверхности в единицу времени.

Стадию химической реакции можно записать в общем виде:

п J(+ 1

ZM ^V IM j=i..r (1.1.1)

Здесь величины а^, и ß^, - это стехиометрические коэффициенты веществ A¿; п- число веществ, г - число элементарных стадий в механизме реакции. Если реакции элементарные, то существуют ограничения на элементы стехиометрической матрицы. Значения otji и ßji могут быть только от 0 до 3. Здесь величины а1, и ß1, - это стехиометрические коэффициенты исходных веществ A¡; и продуктов реакции Bí; соответственно; NA, NB- число исходных веществ и продуктов реакции соответственно.

Если реакции элементарные, то существуют ограничения на эти величины. Значениями a¡ и ßj могут быть только от 0 до 3. Для итоговых уравнений таких ограничений нет.

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

Уравнение стадии гетерогенной каталитической реакции имеет вид :

Ыл Nx Ns NT

Zoe А + >£ ДД + £ /И (1.1.2)

7 = 1 7 = 1 7 = 1 ¡ = 1

где Аг и Bi - соответственно исходные вещества и продукты, находящиеся в газовой фазе; ai, Pi - их стехиометрические коэффициенты; Х{ и Yt-соответственно поверхностные вещества; a'j и p'i - их стехиометрические коэффициенты; NA, Nx ,Nb, Ny - числа соответствующих компонентов реакции. Уравнение (1.1.2) обычно принимает вид:

аА + <—>/» + £ 0•1-3)

i=i i=i

где a, Р - величина, равная 1 или 0, т.е. предполагается, что элементарная каталитическая реакция осуществляется с участием только одной молекулы газообразного вещества, либо без участия этих веществ. Относительно стехиометрических коэффициентов a'i и (3', предполагается, что их величины 1,2,3 (редко) и

Nx Nr

I\a„Zfit*3

i=1 1=1

Скорость стадии вида (1 1 1) и (1.1.2) определяется через разность скоростей прямой и обратной реакций, w+ и w": w=w+-w~ (1.1.4) В равновесии w = 0 и w+ = w'.

Для стадии выполняются соотношения - для гомогенных реакций : w = w+ - w" = -l/(ai*V) * (dNAi/dt) = 1/(P,*V )*(dNBi,/dt) (1.1.5а) где Naj и Nbi - количество вещества системе, моль; V—объём системы, см3;

для гетерогенных каталитических реакций :

+ _ -1 dNAt 1 dNBi

w = w -w =--'- =--L (1.1.50)

a,Sk dt fiSk dt v

где Sk - поверхность катализатора, на котором идёт превращение, см2. Для стадий, идущих без изменений числа молей, соотношения (1.1.5а) и (1.1.56) принимают распространённый вид:

м>

: М> ~М> = ■

■ 1 асА{ = 1 (1св> а, Ж Д ¿//

(1.1.6)

где Сд!, Св1 - концентрации веществ, моль/см3 сек или моль/см2 сек. Согласно закону действующих масс ;

^=к+с:\с%=... = к+1\с2 (1.1.7)

1=1

»-=к-с*с*=... = к-1\с* (1.1.8)

;=1

Коэффициенты к+ и к * - константы скорости реакции, прямой и обратной. Эти константы определяют удельную скорость реакции - скорость при единичных концентрациях реагирующих веществ. Константы экспоненциально растут с увеличением температуры Т.

Стехиометрические коэффициенты веществ, участвующих в реакции сведены в стехиометрическую матрицу Г (з*К), где б - число реакций, N - число реагентов системы. Стехиометрические уравнения сложной реакции могут быть представлены в виде Га = 0, (11.9)

Где а - вектор столбец реагентов. Ограничения, накладываемые на матрицу Г:

Закон постоянства массы атомов по стадиям ГА = 0, (1.1.10) Закон сохранения массы по стадиям ГАМа = ГМ = 0, (1.1.11) Где А(М*т) - молекулярная матрица (т - число химических элементов, входящих в состав реагентов), М - вектор-столбец молекулярных весов веществ, МА - вектор-столбец атомных весов.

Ранг матрицы Г, т.е. максимальное число линейно независимых строк или столбцов, никогда не может быть выше N - т. Это выполняется в силу того, что всегда существует т линейно зависимых столбцов матрицы Г, задаваемой уравнением (1.1.10), гяГОТ-т. (1.1.12)

Ранг матрицы А определяется так:

Rg A=min (N,m). (1.1.13) Как правило rg А = ш. Тогда rg T<N-rgA. (1.1.14)

Это соотношение называется стехиометрическим правилом Гиббса. Вещества, соответствующие линейно независимым столбцам матрицы Г, называются ключевыми, а остальные - неключевыми.

Статическим называется лабораторный каталитический реактор, представляющий собой закрытую систему. Кинетическая модель сложной реакции, осуществляющейся здесь, представляется : С = rTw(c), (1.1.15)

где с - вектор-столбец концентрации веществ, ГТ- транспонированная стехиометрическая матрица; w(c) - вектор-столбец скоростей стадий. Это нестационарная кинетическая модель [ 1 ].

1.2. Методы решения прямых задач для СДУ ХК. Основные трудности Для изолированной многокомпонентной системы, в которой имеет место сложный химический процесс, система уравнений с учетом ОКТ и кинетического закона (1.1.4) может быть записана в виде системы обыкновенных дифференциальных уравнений:

У = ИуЛу{0) =y\t s[0,co) (1.2.1)

Правые части (1.15) являются гладкими (т. е. сколь угодно раз дифференцируемыми) в окрестности начального вектора у°, откуда следует, что решение у(у°, t) начальной задачи Коши существует и непрерывно зависит от координат начальной точки. Это означает, что начальная задача Коши поставлена корректно. Известно также, что решение системы кинетических уравнений (1.2.1) является устойчивым и асимптотически устойчивым по Ляпунову [17,18]. Исходная дифференциальная система как правило, нелинейна, не имеет замкнутого аналитического решения, и для

нахождения решения необходимо применять численные методы. С этой целью исходную систему (1.2.1) аппроксимируют разностной схемой вида

^Ь^ЖО (1.2.2)

И

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

= (1.2.3)

п

где а(Ь) -» 0 при И 0.

Система (1.2.3) является системой преобразованных уравнений с параметром Ь, и очевидно, что (1.2.1) и (1.2.3) эквивалентны лишь в пределе (И -> 0). Однако, в отличие от (1.2.1) нельзя утверждать, что любое решение (1.2.3) вида у(у°,1,11) устойчиво при всяком значении И, Иными словами, процедуры аппроксимации и линеаризации (если последняя проводится) можно рассматривать как источники возмущения на исходную задачу, и, если не принять специальных мер, то решение преобразованной задачи может потерять устойчивость, свойственную решению исходной задачи. Это — первая особенность решения уравнений химической кинетики. Вторая особенность состоит в том, что решения должны иметь свойства положительности и балансности, поскольку ОКТ указывают, что концентрации компонентов ни при каких условиях не могут быть отрицательными и масса системы сохраняется. Проблему устойчивости и положительности решения можно рассмотреть на примере распада одного компонента, т.е. исследуем поведение решения одиночного уравнения

^-ш(1,у)*у,у(0) = у°^(1,у)>0,1е [О,«) (1.2.4) Простейшая аппроксимация уравнения (1.2.4) явной разностной схемой Эйлера первого порядка точности приводит к преобразованному уравнению вида

УП+1=УП[1 - Муп, О ] (1.2.5)

Если шаг 11 столь велик, что хотя бы в одном узле сетки выполняется > 1/Ъ, то численное решение (1.2.5) знакопеременно (или просто отрицательно), что физически бессмысленно. Для одного уравнения эту трудность легко преодолеть, так как за малое время I ~ 1/\у число узлов сетки N = 1 / п сравнительно невелико.

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

у= Ау, \ = 1,2, собственные значения Х[ = 10° и %2 = Ю3- Общее решение,

имеет вид

я

у = ехр(-Л,Л), { = 1,2 [19] . Можно выделить две области решения, в

У=1

которых вклад обоих членов неодинаков. В первой из них, где велик вклад члена с Х\, шаг, определяемый из условия (1.2.4), должен быть равен ХЪ ~ 1,11 ~ 1, а во второй, где велик член ехр( -А^Ь), шаг должен быть ~103 . И, хотя в первой области вклад второго члена незначителен, на всем интервале интегрирования из соображений устойчивости необходимо удерживать шаг именно ~ 103, а не ~1.

В реальных кинетических системах подобная разномасштабность может достигать величин 1015. Физически это означает, что ьй компонент уже почти

достигает своей около равновесной области, в то время как (1 + 1)-й — еще нет. Например, в системе Н21 — О2 концентрации исходных веществ практически не начинают меняться, в то время как концентрации промежуточных веществ уже близки к своим максимальным значениям. Именно такая ситуация наиболее характерна для кинетических систем. С - вычислительной точки зрения это означает, что интервал интегрирования определяется значением самого малого коэффициента скорости, а шаг интегрирования — значением самого большого коэффициента скорости, и

общее число шагов пропорционально их отношению, т. е. алгоритм решения становится недопустимо неэкономичным [20, 21].

Формально система уравнений (1.2.4) квалифицируется как «жесткая», если собственные значения матрицы Якоби 1( у, t) =

|| dfh/dyk II, i, h = 1,...,N, вычисленные на произвольном частном решении, удовлетворяют условиям

Re [Ai(t)] < 0, max {Re [^(t)]} » min { Re [A*(t)]}, t e [0, oo). «Жесткость» системы (1.2.1) на всем интервале решения означает, что постоянно существует ситуация, когда малые изменения значения у\ влекут за собой большие изменения значений fk. В этом случае для правой части системы постоянная Липшица оказывается очень большой по модулю величиной:

|f(y,t)-f(yM|<L|y(t)-y*(t)|.

Малость допустимой величины шага интегрирования h как раз и объясняется тем, что в традиционных методах интегрирования необходимо выполнение условия Lh ~ а (как правило, а ~ 1) с тем, чтобы аппроксимировать как слабо, так и сильно затухающие члены. Итак, при решении прямой кинетической задачи основные трудности вызываются плохой обусловленностью матрицы Якоби. Существуют два принципиальных подхода к решению этой проблемы. Первый из них состоит в том, чтобы «развязать» исходную систему взаимосвязанных уравнений и превратить ее в систему несвязанньк одиночных уравнений, каждое из которых затем решается отдельно. В этом случае не возникает проблемы выбора шага, так как говорить об определении жесткости» для системы (1.2.1) не имеет смысла, и каждое уравнение решается со своим шагом.

Второй подход состоит в том, чтобы тем или иным образом регуляризовать v исходную систему, т. е. преобразовать ее к устойчивой системе, к которой затем можно применить «обычные» численные процедуры (Рунге - Кутта, Адамса и т. д.). На практике такая регуляризация проводится либо на предварительной стадии, и тогда «обычные» численные процедуры

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