Компьютерная обработка химических экспериментов при решении обратных кинетических задач на основе параллельных вычислений Текст научной статьи по специальности «Физика»

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

Вестн. Ом. ун-та. 2012. № 2. С. 194-199.

УДК 519.876.5

М.В. Тихонова, И.М. Губайдуллин

КОМПЬЮТЕРНАЯ ОБРАБОТКА ХИМИЧЕСКИХ ЭКСПЕРИМЕНТОВ ПРИ РЕШЕНИИ ОБРАТНЫХ КИНЕТИЧЕСКИХ ЗАДАЧ

НА ОСНОВЕ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ

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

Ключевые слова: компьютерное моделирование, глобальная минимизация, химическая кинетика, обратная кинетическая задача, кинетические параметры.

Введение

Результат моделирования всегда существенно зависит от адекватности математической модели изучаемому технологическому процессу. При отдельной идентификации математической модели химической реакции для различных температур, при которых проводится натурный эксперимент, возникает проблема учета структурных связей между соответствующими обратными задачами химической кинетики. Для одной и той же реакции при разных температурах, согласно закону Аррениуса [1], константы скоростей элементарных стадий растут с ростом температуры. Один из возможных способов учета таких связей - это компьютерная обработка всей экспериментальной информации путем объединения обратных кинетических задач в одну общую.

1. Математическое описание кинетических задач

Построение кинетических моделей любых химических реакций требует решения двух фундаментальных задач: 1) экспериментальное исследование химического объекта; 2) математическая обработка экспериментальных данных с целью идентификации модели [1].

На этапе математической обработки экспериментального материала производится постановка и решение прямой и обратной кинетических задач.

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

Обратная кинетическая задача - это задача восстановления по экспериментальному материалу вида кинетической модели. Поиск кинети-

© М.В. Тихонова, И.М. Губайдуллин

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

1 Р п Н

Р

к = и = 11 = 1

Р - э ХкИ ХкИ

Ш1П, (1)

где у - вектор минимизируемых параметров; хРы - расчетные значения концентраций наблюдаемых веществ, получаемые в результате решения прямой кинетической задачи на наборе констант у, мольные доли; хЭы - экспериментально полученные значения концентраций наблюдаемых веществ, мольные доли; Н - количество наблюдаемых веществ; п - количество точек эксперимента; Р - количество экспериментов.

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

(2)

У = I к ., к . 1 ] -]

где у - вектор минимизируемых параметров функционала (1); к], к-] - приведенные константы скоростей прямой и обратной ]-й элементарной стадии для одной температуры, 1/ч, соответственно (к, к-] > 0); ],-] - индекс прямой и обратной элементарной стадии.

2. Методы нахождения энергий активаций элементарных стадий

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

Зависимость константы скорости химической реакции к от температуры Т определяется уравнением Аррениуса [1]:

Еа

к = лУнт , (3)

где к - приведенная константа скорости элементарной стадии, 1/ч; Еа - энергия активации, ккал/моль; К - универсальная газовая постоянная, ккал/(моль-К); А - приведенная частота столкновений реагирующих молекул, 1/ч, Т - температура, К.

2.1. Метод наименьших квадратов

Если экспериментальные данные представить в координатах 1п(к) = /(1/Т), тангенс угла наклона полученной прямой линии окажется равным -Еа/К

Е 1

По существующим методикам [3] прямая зависимость 1п(к) = ./(1/Т) строится методом наименьших квадратов (МНК) по имеющимся при различных температурах константам скоростей элементарных стадий, которые находятся путем отдельного решения обратных кинетических задач (1)-(2).

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

2.2. Автоматическое объединение обратных кинетических задач

Один из возможных способов учета структурных связей внутри вычислительной системы - это построение «агрегированных обратных кинетических задач» [4]. Термин «агрегированная обратная кинетическая задача» введен авторами статьи для подчеркивания автоматического объединения всех обратных кинетических задач, необходимых для исследования реакции при различных температурах, в одну общую.

Можно провести непосредственное объединение обратных кинетических задач (1)-

(2), соответствующих различным температурам, в одну:

( т т т т т т \

у = к1, к\, к1.2, А..., ктр, кр.

] -У ] - У ] -]

(5)

1пк = 1пЛ -:

я т

(4)

где у - вектор минимизируемых параметров функционала (1); ],] - индекс прямой и обратной элементарной стадии; Р - количество

химических экспериментов;,7 = 1,Р - температура, при которой проводился г-й хит т

мический эксперимент, К; к _} , к ] - приведенные константы скоростей прямой и обратной ]-й элементарной стадии, соответствующих г-му эксперименту, 1/ч (к, к-] > 0). Но такой способ агрегации обратных кинетических задач влечет за собой увеличение размерности многоэкстремальной задачи

(1), (5) в Р раз. Общая размерность задачи будет равна Р ■ N где N - суммарное количество констант скоростей элементарных стадий, прямых и обратных (размерность одной обратной кинетической задачи для одной температуры).

Во избежание этого, в рамках данной работы глобальный минимум предлагается искать не по константам скоростей, а по параметрам Еа и 1п(А) прямой зависимости Щк) = /(1/Т). Фактор частоты А для простых молекул имеет величину порядка 1013-1015, для сложных - значительно меньшую [1].

Эмпирические значения энергии активации Еа для различных реакций лежат в пределах 0-100,15 ккал/моль, для реакций металлокомплексного катализа - в пределах 10-30 ккал/моль.

В результате для конкретной реакции, не зависимо от количества температур, при которых проводился натурный эксперимент, размерность многоэкстремальной задачи будет оставаться постоянной и равной 2И. Таким образом, семейство отдельных обратных задач (1) преобразуется к одной агрегированной обратной задаче вида:

У =

УІ ’ Уу + *

І = 1, *,

У ■ = Е ,у . ,, = 1пI А .

І а і І + * ^ І

(б)

3. Генетический алгоритм для решения обратных кинетических задач

Решение обратных кинетических задач основано на многократном решении вычислительно трудоёмких прямых кинетических задач с минимизацией функционалов (1), (2) или (1), (6). Это создает основу для эффективного распараллеливания вычислительного процесса.

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

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

ского алгоритма для решения обратных задач химической кинетики [5], включающая в себя 5 стадий (рис. 1).

1. На стадии начального заполнения происходит псевдослучайное заполнение векторов параметров (2) или (6), распределение их между узлами многопроцессорной вычислительной системы (МВС) и расчет функционала невязки (1), (2) или (1), (6) на каждом процессоре.

2. Скрещивание - поиск нового вектора параметров Усгопо образу лучших особей У’ и У":

Усго** = (1 - а)У'+ аУ", 0 < а < 0. (7)

Мутация - отклонение с заданной вероятностью р от полученного в (7) вектора

усгаэв:

Ути = Уст™(1 ± р/100), р - %.

Распределение скрещенных и мутированных точек между процессами МВС и расчет функционала невязки (1), (2) или (1),

(6) на каждом процессоре.

3. Селекция - отбрасывание половины точек с наибольшими значениями функционала.

4. Формирование нового поколения -отклонение на заданный процент от оставшейся половины точек:

г*/2+1

= У = У 41 ± р/100),

у* = у*/2

= У*/2(1 ± р/100).

Распределение новых точек между процессами МВС.

5. Проверка осуществляется по значению минимизируемого критерия: сигналом к остановке процедуры является прекращение уменьшения значения функционала невязки (1).

Рис. 1. Схема генетического алгоритма

На каждой заданной итерации реализовано локальное уточнение полученного минимума методом Хука-Дживса.

4. Сравнение методов определения кинетических параметров

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

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

акция карбоалюминирования олефинов проводилась в присутствии катализаторов Ь^гС12 (Ь = Ср, СрМе5, Ср=С5Н5) при трех температурах (15 °С, 22 °С и 30 °С). В работах [2-3] по экспериментальным данным при различных температурах двумя описанными подходами проводилось моделирование реакции при различных катализаторах. Проведем сравнительный анализ результатов расчета для катализатора Ср^гС12.

В работе [2] проводилась идентификация реакции карбоалюминирования олефинов по отдельности для каждой температуры. Рассчитанные константы скоростей элементарных стадий представлены в табл. 1.

После линейной аппроксимации констант скоростей элементарных стадий методом наименьших квадратов получаются кинетические параметры, представленные в табл. 2.

Для 3-й и 4-й стадий применение отдельного решения обратных кинетических задач для различных температур с последующей линейной аппроксимацией МНК не позволило получить кинетические параметры, поскольку энергии активации попали в отрицательные области (см. табл. 2, рис. 3).

ГГ [Лс]

[Д^

47 (А]Я^ [AJ

Рис. 2. Реакция карбоалюминирования олефинов

Таблица 1

Константы скоростей стадий реакции карбоалюминирования олефинов при Ср22гС!2

№ стадии №, 1/ч

15°С 22°С 30°С

1 прямая 3,65Е+03 9,63Е+02 6,21 Е+03

2 прямая 5,30Е+03 1,22Е+03 5,76Е+03

3 прямая 3,03Е+03 2,22Е+03 2,30Е+00

4 прямая 9,56Е+03 9,64Е+03 7,78Е+03

5 прямая 2,64Е+03 9,06Е+03 6,32Е+03

6 прямая 3,83Е+00 3,52Е+00 5,30Е+03

7 прямая 5,81 Е+03 5,28Е+03 8,27Е+03

8 прямая 4,52Е+02 1,98Е+02 2,74Е+03

1 обратная 2,58Е+03 7,49Е+03 4,86Е+03

2 обратная 7,70Е+03 1,15Е+03 8,60Е+03

Ф= 1,61 1,62 2,52

Для реакций металлокомплексного катализа характерны значения энергий активаций Еа в пределах 2-40 ккал/моль. Помимо того, для 6-й стадии значение энергии активации превышает допустимое, применение МНК приводит к перерасчету полученных при минимизации констант скоростей стадий, что в действительности в 2-3 раза ухудшает полученные минимальные значения целевой функции ф (Н: 6 стадия на рис. 3).

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

решении агрегированной обратной кинетической задачи, приведены в табл. 3.

5. Анализ эффективности параллельного генетического алгоритма

Анализ эффективности использования параллельных алгоритмов обычно состоит в оценке получаемого ускорения Бр и эффективности Ер процесса вычислений [7]. Оценка качества параллельных вычислений предполагает знание наилучших (максимально достижимых) значений показателей ускорения и эффективности, однако получение идеальных величин Бр=р и Ер=1 может быть обеспечено не для всех вычислительно трудоемких задач. Так, достижению максимального ускорения может препятствовать существование в выполняемых вычислениях последовательных расчетов, которые не могут быть распараллелены.

Таблица 2

Кинетические параметры реакции карбоалюминирования олефинов при Ср22гСЬ

№ стадии Ш, 1/ч Еа, ккал/моль А, 1/ч

15°С 22°С 30°С

1 прямая 2,10Е+03 2,77Е+03 3,73Е+03 6,642 2.35Е+08

2 прямая 3,14Е+03 3,33Е+03 3,55Е+03 1,428 3.82Е+04

3 прямая 9,30Е+03 2,85Е+02 6,45Е+00 -83,94 1.45Е-60

4 прямая 9,86Е+03 8,92Е+03 8,01 Е+03 -2,40 147.74

5 прямая 3,49Е+03 5,24Е+03 8,16Е+03 9,808 9.96Е+10

6 прямая 1,07Е+00 3,60Е+01 1,64Е+03 84,650 2.38Е+64

7 прямая 5,29Е+03 6,29Е+03 7,59Е+03 4,154 7.61Е+06

8 прямая 2,49Е+02 6,03Е+02 1,58Е+03 21,341 4.15Е+18

1 обратная 3,35Е+03 4,49Е+03 6,17Е+03 7,066 7.87Е+08

2 обратная 3,91 Е+03 4,23Е+03 4,61 Е+03 1,896 1.08Е+05

ф= 2,94 4,28 4,07

фср = 3,76

1/Т 1/Т

Рис. 3. Прямая зависимость 1п(к) = /(1/Г) для 3-й (слева) и 6-й (справа) стадий

Таблица 3

Кинетические параметры элементарных стадий схемы 1 реакции карбоалюминирования олефинов

в присутствии катализатора Ср22гСЬ

№ стадии Ш, 1/ч Еа, ккал/моль А, 1/ч

15°С 22°С 30°С

1 прямая 1,91Е-02 3,30Е-02 5,96Е-02 13,06 1,72Е+08

2 прямая 3,40Е+02 4,94Е+02 7,41 Е+02 8,97 2,33Е+09

3 прямая 1,78Е+03 2,51 Е+03 3,65Е+03 8,27 3,57Е+09

4 прямая 1,88Е+10 3,14Е+10 5,48Е+10 12,33 4,69Е+19

5 прямая 4,28Е+06 5,54Е+06 7,33Е+06 6,18 2,20Е+11

6 прямая 7,43Е+05 8,37Е+05 9,53Е+05 2,86 1,11Е+08

7 прямая 5,49Е+07 6,06Е+07 6,75Е+07 2,38 3,57Е+09

8 прямая 5,91 Е+03 7,40Е+03 9,46Е+03 5,42 7,99Е+07

1 обратная 3,20Е+01 4,32Е+01 6,00Е+01 7,21 9,94Е+06

2 обратная 2,60Е+07 3,06Е+07 3,65Е+07 3,94 2,63Е+10

2,67 2,46 3,18

фср= 2,77

Рис. 4. Анализ эффективности генетического алгоритма

На кластере Башкирского государственного университета произведено тестирование реализованного параллельного генетического алгоритма для определения кинетических параметров реакции карбоалюми-нирования олефинов (обратная кинетическая задача (1),(2) или (1),(6)) в присутствии катализатора (CpMe5)2ZrCl2 при 15 °С [5], исследуемой в Институте нефтехимии и катализа РАН. Для количества 1-28 процессоров при начальной популяции 400 и популяции скрещивания 28 особей произведено измерение времени выполнения одной итерации алгоритма. Рассчитаны получаемое ускорение и эффективность выполнения программы (рис. 4).

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

Заключение

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

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

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

Внедрение и тестирование вышеуказанных методов проводилось для определения кинетических параметров новых реакций карбоалюминирования и циклоалюми-нирования олефинов [7]. Новые методы позволили существенно сократить время построения кинетических моделей.

ЛИТЕРАТУРА

[1] Царева З. М., Орлова Е. И. Теоретические основы химической технологии. Киев : Вища школа. Головное издательство, 1986. 271 с.

[2] Губайдуллин И. М., Рябов В. В., Тихонова М. В. Применение индексного метода глобальной оптимизации при решении обратных задач химической кинетики // Вычислительные методы и программирование. 2011. Т. 12. С. 137-145.

[3] Спивак С. И., Губайдуллин И. М., Вайман Е. В. Учебное пособие. Уфа : Редакционно-издательский отдел Башкирского государственного университета, 2003. 110 с.

[4] Тихонова М. В., Рябов В. В. Решение агрегированных обратных задач химической кинетики параллельным индексным методом глобальной оптимизации // Научный сервис в сети Интернет: экзафлопсное будущее : труды Международной суперкомпьютерной конференции. М. : Изд-во МГУ, 2011. С. 572-579.

[5] Тихонова М. В., Губайдуллин И. М., Лаврентьева Ю. С., Масков Д. Ф. Распараллеливание агрегированных обратных кинетических задач математического моделирования реакций металлокомплексного катализа // Системы управления и информационные технологии. 2011. № 4 (46). С. 10-14.

[6] Parfenova L. V., Gabdrakhmanov V. Z, Khali-lov L. M., Dzhemilev U. M. On study of chemose-lectivity of reaction of trialkylalanes with alkenes, catalyzed with Zr n-complexes // J. Organomet. Chem. 2009. V. 694. №. 23. P. 3725-3731.

[7] Гергель В. П. Теория и практика параллельных вычислений. М., 2007. 423 с.