CC BY
20
УДК 533.7/9+541.124]:519.6
ГЕНЕРАТОР ПОДРОБНЫХ КОМПЬЮТЕРНЫХ МОДЕЛЕЙ ДИНАМИКИ РЕАГИРУЮЩЕГО ГАЗА: НОВОЕ ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ СИСТЕМНЫХ ИССЛЕДОВАНИЙ В ТЕХНИКЕ И ЭКОЛОГИИ
© 2006 г. О.В. Яценко, Е.Н. Ладоша
Most of required modifications are made to improve the existent model's generator aimed at computer simulations of flows with chemical reactions. These are concerned with the thermo-chemical self-consistence of elementary acta, with the inhomohenity of reacting gas-flow, with the structural and/or dynamics stability of medium, as well as with the reaction and transport arrangement according to it's general importance. Furthermore, respective empirical database on micro-kinetics is developed essentially. All together improvements made are a perfect tool for numerous detailed researches in modern technology and actual ecology.
Комплексное в физико-химическом отношении содержание задач технической экологии и безопасности жизнедеятельности побуждает использовать разнообразные средства автоматизации вычислительного эксперимента (ВЭ). Поскольку элементами данной предметной области являются перенос вещества, лучистой энергии и химические реакции, понятна острая востребованность специализированных вычислительных пакетов газовой динамики (ГД), химической кинетики (ХК) и радиационного переноса (РП) [1, 2].
История подобных средств автоматизации началась два десятилетия назад, когда «ручное» программирование ЭВМ стало главным фактором, ограничивающим подробность моделей и реалистичность имитаций. Выход из сложившейся ситуации был найден в массовой разработке всевозможных средств автоматизации программирования, предназначенных для интенсификации прикладных исследований. Следующим шагом в развитии научных информационных технологий стала интеграция отдельных средств эффективного программирования - в проблемно-ориентированные информационные комплексы. Главным достижением здесь явилось создание генераторов моделей (ГМ). Этим термином обозначаются программные средства, позволяющие: 1) получать, хранить, сортировать, изменять, уточнять и т.д. теоретические и/или эмпирические данные; 2) агрегировать фактическую, теоретическую и эвристическую информацию, формируя «на выходе» исполняемый компьютерный код; 3) планировать и реализовывать предметный ВЭ.
Проблемы, стоящие перед разработчиками востребованных в технической экологии автоматизированных систем ГД, сводятся к необходимости строить адаптивные сетки - т. е. регулярно осуществлять ре-дискретизацию расчетной области в процессе эволюции течения. Идея такого подхода состоит в локальном сгущении сетки только вблизи развивающихся особенностей: при этом тонкая структура течения разрешается дискретной моделью, а суммарный объем вычислений ощутимо сокращается. Успешной реализации этой идеи, однако, препятствует очень высокая «жесткость» реальных ГД задач: даже локальное измельчение расчетной сетки сопровождается катастрофическим увеличением необходимого объема вычислений.
Первыми в приложениях, включая техническую экологию, оказались востребованными простейшие ИММ химических процессов - модели химически рав-
новесных сред. В них химический состав однозначно определяется содержанием отдельных элементов Q и термодинамическими параметрами - температурой T и давлением P среды. Отсюда актуальность автоматизированных систем (АС) термохимических расчетов, наиболее удачными из которых являются отечественные ИВТАНТЕРМО, TETRAN, АСНИР ТОТ-МГТУ и зарубежные STANJAN и CHEMKIN. Однако пренебрежение динамикой или кинетикой реакций делает равновесное описание недопустимо грубым в современных приложениях, где требуется устанавливать оптимальные соотношения между локализованной нелинейностью и геометрией реагирующей системы. Потребность адекватно описывать реакции, протекающие с конечными скоростями, удовлетворяется моделями химической кинетики, реализованными в форме АС. Первенство здесь принадлежит американцам [3,4]. Однако если учесть затрачиваемые на проектирование и реализацию подобных АС усилия, а также их последующую окупаемость, отечественные системы АВО-ГАДРО [5], БАНКОН [6] и КИНКАТ [7] практически не уступают программам KIVAII [3] и MACRON [4].
Автоматизированные системы ХК [3-7] принципиально позволяют рассчитывать химические равновесия как предельные (в бесконечно отдаленный момент времени) состояния материально и энергетически консервативных реагирующих систем. В этом случае говорят о расчете параметров предельного состояния методом установления. Расчет равновесного состава методом установления требует определенной аккуратности - с целью исключить многозначность модели химической динамической системы (ДС) в результате игнорирования термохимической потен-цильности (материально замкнутых) реагирующих сред. В случае открытых реагирующих систем корректность модели в части потенциальности превращений особенно важна, поскольку является непременным условием корректной идентификации причин и механизмов сложной физико-химической динамики.
Предельные на данном этапе номенклатурные показатели ИММ, составленные при помощи средств автоматизации компьютерного моделирования (КМ), приведены в табл. 1. Как следует из представленных данных, сферой приложения АС в технической экологии и безопасности жизнедеятельности являются ИММ высокого разрешения как по числу существенных переменных, так и в пространственно-временном отношении.
Таблица 1
Современные номенклатурные пределы моделей, описывающих физико-химические процессы в технической экологии и безопасности жизнедеятельности
Процесс Приложение (цели экологии и безопасности) Характеристика ИММ
Течение среды Оптимизация взаимодействия летательных аппаратов с земной атмосферой (снижение энергетических потерь, экономия топлива, повышение надежности) Двумерное (2Б) поле термогазодинамических параметров газа на неравномерной треугольной сетке, состоящей из ~104 узлов
Химические превращения Оптимизация способов и режимов сжигания углеводородных топлив в энергетических установках (повышение КПД, ресурса двигателя, снижение токсичности выхлопных газов) Ящичная модель химического реактора для ~ 104 реакций между ~ 3-103 компонентами
Перенос лучистой энергии Прогнозирование уровня радиационной нагрузки на поверхность спускаемых космических аппаратов (оценка и минимизация рисков, связанных с термическим разрушением конструкций) Одномерное уравнение переноса излучения, осредненное по частоте, температуре излучающего слоя; при осреднении оптических коэффициентов учтены как непрерывны спектр, так и ~ 106 спектральных линий
Основой АС моделирования реакций, реализованной авторами, послужила система КИНКАТ [7]. АС КИНКАТ реализована на 1ВМ-совместимых персональных компьютерах и предназначена для интеллектуализации теоретических исследований физико-химических процессов в газах и плазме. В ее состав входят база данных (БД), система управления базой данных (СУБД) и генератор кинетических схем и машинного кода (ГКСМК). Структурно-функциональная схема автоматизированной системы, а также содержание базы данных и фрагмент пользовательского интерфейса приведены на рис. 1.
Авторы расширили функциональные возможности АС КИНКАТ [7] и восполнили «пробелы» в химическом сегменте БД, выполнили необходимое согласование кинетических и термодинамических коэффициентов. Целями названных действий ставились: 1) «уплотнить» номенклатуру в субмоделях реакций [6,7]; 2) согласовать химический модуль с существующим термохимическим потенциалом; 3) обеспечить возможность автоматизированной генерации кинетических схем реакций в распределенных системах; 4) реализовать автоматизированный контроль устойчивости кинетической схемы к значениям констант скоростей реакций и начальных концентраций реагентов; 5) усовершенствовать систему алгоритмов трассировки эволюционных уравнений, направленных на выявление ведущих процессов.
Накопление, систематизация и согласование кинетических данных. Все реакции рассматриваются как пары взаимно обратных процессов, константы скоростей которых Ь/ и k^° удовлетворяют зависимостям
lnkkf = ln4/ - E//T ,
lnkkb = ln4kb - Ekb/T , (1)
В формуле (3) индексом обозначены стандартные условия - Р51 = 1 атм, Т51 = 298 К. Входящие в (3) величины ЛХ^ ь, ЛСр k и ЛН51 ь удобно представить в векторной форме
ASsfk ^ 's, ^
ACpk = Avk CPki , (4)
AHstk у U у
где .4/, Аьь, Е/ и ЕЬ° - коэффициенты в уравнении Ар-рениуса; энергия активации Еь определена в единицах Я; Ь = 1, . . . М/2 - номер реакции. Отношение констант скоростей в уравнениях (1) равно надежно определяемой экспериментально константе равновесия К
1пКь = (1пА/ - 1п4ьь) - (Е/ - Еь°)/Т , (2)
теоретически корректным шаблоном которой служит 1пКь = (Л^ ь - ЛСр ь) + ЛСр ь 1п(Т/Тй) --(ЛН ь/Т - ЛСр т/Т) . (3)
где Л vjk = Vюь ^'ь - изменение числа атомов /'-го
вещества в Ь-й элементарной реакции; Х/, Ср / и Н/ -соответственно энтропия, теплоемкость и энтальпия участников реакции (в стандартных условиях).
Поскольку величины Х/, Ср / и Н/ надежно измеряются и затабулированы для всех практически значимых веществ, определение неизвестных кинетических коэффициентов Аь и Еь в одном из уравнений (1) по заданным коэффициентам второго уравнения всегда выполнимо на основании термохимических параметров реагентов и продуктов Х/, Ср / и Н/ и уравнений связи (2)-(4). Алгоритм сильно упрощается, если теплоемкости всех компонентов Ср / положить независящими от температуры. Описанное «уплотнение» БД осуществлено для набора реакций [6-8], дополненного сведениями о нескольких десятках элементарных процессов - деструкции высших и реформинга низших углеводородов. Результатом совершенствования БД [6-8], содержащей первичные кинетические данные 0 ~ 600 элементарных процессах, является двукратное расширение номенклатуры учтенных реакций.
Чтобы обеспечить термохимическое согласование реакций в генерируемых кинетических схемах, реализовано следующее алгоритмическое усовершенствование АС КИНКАТ. При наличии в кинетической БД первичных сведений о прямом и обратном процессе составитель-исследователь должен указать, какие из них достоверней: кинетические коэффициенты о другом парном процессе автоматически рассчитываются на основании известных Х/, Ср / и Hi по формулам (2)-(4).
Если же для некоторых реагентов (продуктов) не известны значения Х/, Ср / и Н/, но зато известны кинетические характеристики множества элементарных
ПОЛЬЗОВАТЕЛЬ (инженер-исследователь)
t
интерактивный обме« даннь/ми
t
Головной модуль системы КИНКАТ
Автоматизированная снстена хиническон кннетнкн КИНКАТ
пГ
I Описание возможностей системы Управление данными Генерация кода Выход
дашм I за а Н _§|| о в
"Л _
изучение, возможностей, функций и инт ерфейса систе мы KtftHtA Т
Справочная подсистема:
стр уктура,
управление,
данные,
для псе к учтенном веществ и снами овйк оо сто !ний
База
термодинамических данных
Теплота образования вещества
Частоты молекулярных колебаний a>j
Г
интерактивный обмен данными
Система управления базой кинетических данных
Генератор кинетической схемы
* _ л А
обмен функциональной и содержательной информацией
База исходных кинетических
данных для химических и ф отохимических реакций
й ля учме нных реакций
Модуль факторного анализа
Энтропия образования вещества (7)
Теплоемкость вещества
сАл
База кинетических данных
N.O.NO.N/. О/. N+, 0+, NO+,
Параметрически возмущенный воздух
СН. CN. ОН. HNO .С О,HCl, SH.
Химически возмущенный воздух
Число учтенная, реакиуй:
■ 100
|Ы?М). CO?iDÜ1 V 0?(э). Oi'DV
Фотохимия и перенос квантовых возбуждений hca.....
■ 800
■200
Число учте жых химически*. элеменл)oe -S, веществ и ихвабужденнок состояний ~100
Рис. 1. Структурно-функциональная схема усовершенствованной АС КИНКАТ [7]
реакции с их участием, авторами реализован и включен в состав АС КИНКАТ алгоритм, позволяющий определить согласованный термохимический потенциал для заданной системы химических реакций. Предполагается, что теплоемкости Cpi всех компонентов равны, не зависят от температуры Т, и в каждом элементарном акте суммарная теплоемкость реагентов и продуктов совпадают. Эти предположения не влияют на адекватность большинства кинетических моделей в инженерно-экологических исследованиях и позволяют свести (3), (4) к
\пК = Д^ к- ЛЯ84 к/Т , (5)
ÍДS^ ^
'stk AH stk
= Дк
ik
(st Л
v Hi у
(6)
Удобно предварительно симметризовать задачу, представив температурный ход всех реакций моделируемых реакций в аррениусовском виде (1). Для этого заданные в наиболее общей форме температурные за-
висимости констант скоростей
kk (T) = AkTBke-Ek/T
преобразуются в kk (t)= Ake Ek /T . Естественным условием близости {Ak, Ek} ~ A, Bk, Ek } служит
I
Ake~Ek /T / AkTBke~Ek /T
2
dT ^ min,
(7)
где Ттш -г- Ттах - температурный диапазон «рабочих условий» модели.
Считая неизвестными все Si и Ии определим их методом наименьших квадратов. В качестве целевой функции (ЦФ) выберем минимальную суммарную
T
T
погрешность отклонения синтезируемого термохимического эффекта реакций (ДvikSi, Ду1кЯ} от заданного посредством набора известных констант скоростей реакций (ASst к, ДЯ* k}:
ЦФs = ЕЕ (Д VkS, - Д^ к ) "> min,
к i
ЦФя = ЕЕ(кЯг -ДЯstk) ->min,
(8)
к i
где ASst к = 1о4/ - 1пАД ДЯst к = Е/ - ЕкЬ. Дифференцированием по искомым параметрам Sj и Hj из (8) следует
(9)
El ЕДVkДу}к \Si = ЕДvЗкASstк;
к 1 i У к
Е\ЕДугк Ду}к \Я, =ЕД Vjk ДЯ st к
Обозначая ау = Е Л vik Л у , ¿у = 2 Л vik ь, И = Е Л vik ЛЯ ь, имеем
ах = ж , ах Ил = к , (10)
откуда
= а-1 х ж, Н^ = а" 1 х к . (11)
Условие разрешимости (10) detа Ф 0 не является обременительным и выполняется для большинства «плотных» многокомпонентных систем реакций.
Во многих случаях константы скоростей элементарных реакций экспериментально измерены при некоторой фиксированной температуре (как правило, при 298 или 300 К). В то же время для практического моделирования необходимо знать также их зависимость от температуры. Расчет температурных факторов осуществляется в усовершенствованной системе КИНКАТ, исходя из того, что известное значение константы скорости представляет собой произведение частоты бинарных столкновений, стерического и барьерного факторов. Полагая частоту двойных соударений равной 210-10 см3/с, а стерический множитель ~ 0,1, получаем критическую или максимальную величину константы скорости бимолекулярной реакции ькрит и 210-11 см3/с (в нормальных условиях). Отличие фактического значения от критического относится на счет барьерного или активационного фактора. Его величина определяется выражением типа е ЕА/Ят, а температурная зависимость - выражением
ь(Т) =
где индексом «табл» помечены табличные значения константы и температуры. Очевидно, в такой же корректировке нуждается и температурная зависимость константы обратной реакции: соответствующую величину энергии активации здесь следует увеличить на
-Я Ттабл1п(ь^табл/ь^крит).
Очень важно также, что некоторые элементарные реакции «меняют порядок» при увеличении давления. Конечная скорость внутримолекулярных превращений ограничивает темп соответствующих элементарных реакций при достаточно интенсивном термическом возбуждении реагентов. В результате интенсификация молекулярных столкновений сопровождается равноценным ускорением реакции лишь до определенных пределов по давлению (строго говоря, концентрации активатора), после чего наблюдается понижение порядка реакции со второго - на первый (и с
третьего - на второй для обратных реакций). При генерации моделей реакций усовершенствованной АС КИНКАТ константа скорости таких процессов вычисляется по формуле
ь([М],7) =
=(ь^+1([М],7)[М] ь„(Т))/(ьп+1([М],7)[М]+ь„(Т)), (13) где [М] - концентрация активатора реакции; ь„+1([М],Т) и ^„(Т) - константы скоростей реакций, лимитируемых концентрацией активатора и темпом структурной трансформации активированного комплекса соответственно. Для таких процессов сперва вычисляется константа скорости прямого процесса переменного порядка согласно (13), а затем константа скорости обратного процесса переменного порядка ьь([М],Т), равная К(Т)/ьг([М],Т). Собранные в усовершенствованной БД константы скоростей обратных реакций переменного порядка представляют собой предельные значения при бесконечно малых и, наоборот, при очень высоких значениях [М].
Описанный аспект проектирования моделей химизма тестировался путем сравнения рассчитанных по (13) констант обратных процессов с известными кинетическими данными. Установлено хорошее согласие данных, получаемых описанным и независимыми способами (рис. 2), что свидетельствует о надежности разработанных методик и их программной реализации.
Рис. 2. Сравнение рассчитанной при помощи АС КИНКАТ концентрационной зависимости константы скорости реакции, обратной процессу: а - СН3+Н ^ СН4, с данными [9]; б - СН2О+Я ^ НСО+ЯН с данными [10]. Линиями показаны расчетные, точками - экспериментальные данные. Чтобы повысить наглядность, константа скорости реакции для Я = Н увеличена, а для Я = С2Н5 - уменьшена в 1000 раз (Я - алкильный радикал)
Неизвестные термодинамические характеристики реагирующих веществ, как и кинетические, можно оценить, предполагая их аддитивность с учетом (структурного) подобия химических соединений. «Термодинамический» вектор {Н, Б, СР} не представленного в справочнике [11] вещества (в качестве примера рассмотрим изученное вещество С5Н12) составим как линейную комбинацию известных векторов одного из гомологов этого вещества и пары отличающих эти вещества радикалов. Здесь возможны варианты: от выбора базовых элементов зависит точность оценки. Очевидно, надежность прогноза определяется его «дальностью», т.е. сходством качественных и количественных показателей базового и прогнозируемого соединений, а также выбранных радикалов замещения (см. пример в табл. 2).
Данные на рис. 3 иллюстрируют хорошее воспроизводство констант равновесия важных для приложений реакций с помощью усовершенствованной АС КИНКАТ. Таким образом, результаты всестороннего тестирования усовершенствованной АС свидетельствуют о надежности инструментальных расчетных методик, первичных термохимических и кинетических данных, а также высоком качестве программной реализации.
Таблица2
Техника восстановления термодинамической информации методом подобия
Вариант Способ вычисления H, кДж/моль S, Дж/моль,К с, Дж/моль,К
1 С5Н12 = С2Н6 + +3(СзН7 - С2Н5) -144,1 325,2 113,0
2 С5Н12 = С2Н6 + +ЗХС2Н5 - СНз) -178,5 289,0 96,5
3 С5Н12 = СзН8 + +2"(СзН7 - С2Н5) -143,5 333,7 113,7
4 С5Н12 = С3Н9 + +2ХС2Н5 - СНз) -166,4 309,5 107,7
5 С5Н12 = С4Н10 + +(СзН7 - С2Н5) -144,5 341,9 117,9
6 С5Н12 = С4Н10 + +(С2Н5 - СНз) -156,6 329,8 112,4
Табличные значения [11]: н-пентан 2-метилбутан 2,2-диметилпропан -146,4 -154,5 -166 348,4 343,0 306,4 122,6 120,6 121,6
Скорость и принимает одно и то же значение для всех ячеек в одномерном случае и подчиняется уравнению непрерывности для каждого химического компонента и температуры в при цилиндрической и сферической геометрии задач. Схеме дискретизации (рис. 4) соответствует следующий вид обменных слагаемых иЧвв уравнении (15)
■Ч-1Г/Л в>1
иУд
дискретизация
j-X - urf~ldJ
z z 7 - 7-1
(16)
и большая система обыкновенных дифференциальных уравнений (ОДУ)
dd/dt = f(d) - Z-
'J-1' j
z-w-1
irJ-i
- uJrj ~le]
Рис. 3. Константа равновесия реакции СО + Н2О ^ СО2 + Н2 , рассчитанная при помощи АС КИНКАТ (линия) в сопоставлении с данными [12] (точки)
Генерация моделей распределенных химических реакций. В рамках усовершенствования авторами реализованы два варианта генерации кинетических схем в неточечных системах: 1) модель «проточного реактора» и 2) модель «кинетика-диффузия». В модели проточного реактора фиксируются границы ячеек, на которых задается скорость течения среды. При этом к реакционному слагаемому _Дв) в кинетическом уравнении
йвШ = /в) (14)
генератор кода автоматически добавляет конвективное слагаемое иЧв. В уравнении (14) и далее будем использовать обозначение в= {сь с2, ..., Т}, подразумевая покомпонентное применение оператора V в записи
дв/д/ = /(в) - и Vв /покомпонентно/, (15)
где и - скорость течения среды. Имеется возможность выбрать пространственную размерность задачи С, от 1 до 3. В пределах]-й ячейки газ считается химически и термодинамически однородным и определяется параметрами в == {в^}. Схема дискретизации задачи по пространству показана на рис. 4 а.
zz rj - rj-i
(17)
От (14) система (17) отличается в двух отношениях: во-первых, размерность «химического вектора» в увеличилась в 3 раз, где 3 - число пространственных ячеек; во-вторых, система (17) при нестационарном поле скоростей и = и(/) утрачивает автономность. Оба названных усложнения модели заметно отягощают численное интегрирование (17).
Сведение уравнений с частными производными (УЧП) к большой системе ОДУ методом прямых здесь выполнено в соответствии с формализмом Эйлера: наблюдение за эволюцией среды осуществляется из фиксированных точек пространства (ячеек). Если же реакции протекают на фоне диффузионного переноса, более удобным оказывается формализм Ла-гранжа, состоящий в наблюдении эволюции фиксированных масс (объемов) жидкости или газа. Это обстоятельство учтено при реализации стандартных ситуаций «кинетика - диффузия» усовершенствованным генератором кинетических схем.
В условиях диффузионного переноса реагентов уравнение (14) дополняется слагаемым (Ув
/покомпонентно/)}. Диффузия, в том числе турбулентная, обладает свойством автомодельности или самоподобия: при «высокой» симметрии задачи можно перейти к новым переменным, в которых решение задачи инвариантно. Другими словами, задав некоторое разбиение пространства с химически инертным
диффундирующим веществом (примесью) на элементарные ячейки Т/, можно определить закон разбухания этих ячеек т3(/), при котором количество примеси в каждой из ячеек сохраняется сколь угодно долгое время. Если же примесь вступает в химические реакции и ее концентрация меняется, между соседними ячейками возникают дополнительные (градиентные) потоки вещества. Схема дискретизации задачи «реакция - диффузия» в усовершенствованной АС КИН-КАТ показана на рис. 4 б.
После дискретизации УЧП «кинетика - диффузия» согласно рис. 4 б
дв/дГ = /(в) - <у{£(Ув) /покомпонентно/} , (18) преобразуется в большую систему ОДУ
в _в3-1 )Т/_-1(()
+ Z
dQ/dt = f(93) + f
j -Qi yf-(()
f f ri - rj-1
f f rf - rf-1
'i-i dt
dr
i dt
(19)
Непосредственная подстановка показывает, что дискретная модель (19) материально консервативна: Лив 3= - д\яУ3, где V3 - объем/-й ячейки.
Оценим влияние усложнения модели реакций, обусловленное «наложением» транспортных эффектов, на трудности численного интегрирования соответствующих динамических уравнений. Как и предшествующие, усовершенствованная версия АС КИНКАТ использует интегратор больших систем нелинейных ОДУ абсолютно устойчивым методом Гира [13].
О1
. = 1
О-/0 rj(0
изотермических газах температура рассматривается как отдельный компонент в числе I).
Успешность интегрирования ограничивается, главным образом, двумя факторами: «жесткостью» системы и разреженностью якобиана. Оба эти фактора усугубляются при увеличении числа различаемых реагентов (и продуктов) I, а также при повышении пространственно разрешения модели (числа ячеек I). Сравнить структуру якобиана точечной и соответствующих распределенных систем с конвекцией и диффузией позволяет схема рис. 5.
Нетрудно подсчитать, что общее число динамических переменных в уравнениях (17), (19) равно I х I, доля ненулевых элементов д/(в)/дв в каждой пространственной ячейке и I-12, а в совокупной матрице и 2 -г 3 Г112!1. При типичной жесткости химических систем в десять порядков величины, которая увеличивается при учете пространственной неоднородности реагирующих (течений) газов, возможности инструментального интегратора ограничиваются кинетическими схемами размерности/п~3< 100.
4
и
■ ■ ■
1 ■
1/г ■
"ХИМ ■
i - номер вещества (i = 1, ... I)
И
]/г
V
1/ги
'misiI
е
Vr 1 мм
•/г:
Рис. 4. Разбиение реагирующего газа на ячейки в предусмотренных усовершенствованной АС КИНТАТ стандартных моделях с конвективным (а) и диффузионным (б) переносом
Ключевую роль в устойчивых (неявных) методах интегрирования ОДУ играет техника обращения якобиана системы д/(в)/дв. Данная процедура требует порядка 13 операций, для матрицы д/(в)/дв размером I х I отвечающей точечной системе I веществ (в не-
м'"Н>|
Тш'И'П
Ш
(л;и|и[)|
'/г.
ншая
Li.
1/г
D
) номер пространственной ячейки (/ I. ... ,/) I /*/+<- сквозной номер «локализованного вещества»
Рис. 5. Структура «химического» якобиана: а - в точечной постановке; б - в системе с конвективным переносом; в - в системе с диффузией
а
а
б
б
в
Это ограничение не препятствует успешному моделированию на ЭВМ многих (но отнюдь не всех) актуальных задач технической экологии. Отметим, что погрешность аппроксимации реалистичных УЧП системами ОДУ (17), (19) имеет первый порядок точности по пространству (при равномерном по массе разбиении (19) - второй) и переменный (выбирается интегратором [13] от 1 до 5 в зависимости от гладкости якобиана) по времени.
Следует отметить особо, что автоматическое включение в кинетическую схему обратных реакций обеспечивает непременное заполнение главной диагонали д/(9)/дв1 якобиана: в результате улучшается его обусловленность и повышается устойчивость машинного интегрирования ОДУ.
Исследование структурной и динамической устойчивости моделей. Принципиально важной при моделировании течений реагирующих газов является структурная устойчивость совокупных моделей и их отдельных блоков. Даже точечные кинетические модели представляют собой сильно нелинейные ДС с очень широким частотным спектром. Материально и энергетически замкнутые параметрически однородные объемы реагирующих газов в динамике диссипа-тивны: аттракторами служат равновесные состояния. Подобные системы, безусловно, являются грубыми при разумном «шевелении» их внутренних констант (скоростей реакций).
Усложняя постановку и приближая условия реакций к действительности, приходится в различной степени «открывать» и рассредоточивать модельную систему ХК. Вопросы структурной устойчивости таких реалистичных ИММ остаются открытыми. Некоторыми гарантиями здесь можно заручиться, если воспользоваться следующим приемом. Задача решается дважды: сперва с базовым комплектом констант реакций {к/, кД к/, к2ь, ... }, а затем - с возмущенным {к/, кД к/, к2ь, ... }дзЬ в котором каждое значение к■ умножается (а к,ь соответственно делится) на случайное число, логарифм которого равномерно распределен в интервале 1/2 -г- 2 или 1/3 -г- 3. Полученные решения сравниваются на предмет взаимного сходства (в части целей моделирования). Если различие базового и возмущенного решений соразмерно с уровнем возмущения, модель признается проблемно или условно устойчивой. Интересно при этом также выявить компоненты и процессы в кинетической схеме, наиболее чувствительные к возмущению констант. В проблемно устойчивых моделях легко возмущаемые описанным способом компоненты не важны для совокупного химизма: их можно исключить из рассмотрения.
Решению проблемы структурной устойчивости кинетических уравнений служит отсутствие математически безупречных методов численного интегрирования жестких нелинейных ОДУ (для УЧП эта проблема столь сложна, что здесь даже не обсуждается). Наиболее употребительный при моделировании реакций интегратор Гира [13] работает по принципу сужения частотного спектра: постепенно из числа учитываемых в область «машинного шума» выбрасываются все более медленные процессы, что позволяет осуществлять прогрессирующий «разгон счета». Эффективная для замкнутых точечных реакторов техни-
ка [13] не срабатывает в случае недостаточно гладких и тем более разрывных потоковых слагаемых в кинетических уравнениях. Хотя в этом случае проблема разрешима - через перезапуск интегрирования в точках (в моменты) нарушения гладкости внешних потоков. Подчеркнем, что интегрирование систем с обостряющимся химическим ядром якобиана, в которых принципиально возможен режим «жесткой турбулентности» и/или «переключающей перемежаемости», проблематично даже с помощью такого робаст-ного алгоритма, как [13].
Есть и другой способ возмутить внутренние параметры кинетической ДС и определить при этом не только структурную устойчивость модельных уравнений вообще, но также их устойчивость в условиях конкретной неопределенности наличных констант скоростей реакций. Его реализовать особенно просто, если в распоряжении исследователя имеются независимые экспериментальные значения констант для всех прямых и обратных реакций. Также как и в случае стохастического возмущения, необходимо выполнить два сравнительных интегрирования модельной ДС химизма: в первом кинетическая схема снабжается известными константами «прямых» реакций и термодинамически согласованными (рассчитанными) константами реакций «обратных», во втором случае поступают наоборот. Сравнение полученных результатов позволяет судить о практической устойчивости ИММ.
Более последовательным способом установить грубость кинетических моделей является реализация ансамбля специальным образом возмущенных имитаций с общими начальными условиями и построение на их основе доверительных (вероятностных) коридоров для каждой динамической переменной. Определяемые так коридоры не имеют четких границ в отличие от устанавливаемых путем сравнения только двух «граничных» случаев - невозмущенного и возмущенного (рис. 6).
_ а
граничные кривые
М
Рис. 6. Оценка достоверности кинетической кривой методами «граничных» (а) и «стохастических» (б) возмущений. Символом Р(с) обозначена плотность вероятности того, что концентрация г'-го вещества равна сг
Когда приемлемая структурная устойчивость ХК-блока установлена, необходимо оценить влияние начальных условий - химических и термодинамических на результат моделирования. Реализованное авторами усовершенствование АС КИНКАТ позволяет автоматизировать случайное возмущение начальных условий (НУ) в заданных границах и установить пределы прогноза эволюционной динамики для компонентов в/ со сложной динамикой.
Лишь после выполнения описанных здесь верификационных процедур можно судить об адекватности модели, ее надежности, диапазоне работоспособности. Новые средства генерации и анализа моделей в составе АС КИНКАТ создают условия, когда рутинные исследования динамических качеств модели (в частности, течений газа с реакциями) перестают быть непреодолимым препятствием в получении новых теоретических результатов.
Выделение ведущих химических реакций в большой совокупности. Располагая детальной моделью реагирующего неоднородного газа, можно решать задачу построения минимальной модели превращений. Минимальная модель реакций включает наиболее «быстрые» или далекие от вырождения концентрационные переменные. Хотя современные теоретические заделы, экспериментальные данные и информационные технологии позволяют конструировать очень громоздкие кинетические схемы превращений, их пока не удается использовать в качестве блоков интегрированных проблемных моделей (например, в сочетании с одно- и тем более 2 -г 3-мерной моделью движения среды).
Поэтому весьма привлекательно заблаговременно выделить из всей совокупности реакций небольшое подмножество, отвечающее за «важные черты» подробной кинетики. Составляющие усеченный блок реакций элементарные процессы назовем определяющими, если усечение схемы превращений проведено в характерных условиях основной задачи. Так получаемые минимальные модели химизма без особых проблем сочетаются с аналогичными (асимпотическими) моделями прочих эффектов в больших междисциплинарных ИММ. Описанный прием позволяет успешно решать сложные задачи динамики реагирующих газов, в том числе технические и экологические.
Методы важностного ранжирования элементарных реакций можно разделить на две категории: 1) анализа чувствительности и 2) целевого поведения. В основе анализа чувствительности лежит гипотеза, согласно которой изменение совокупной динамики системы в сильной мере определяется изменением важных элементов и слабо зависит от элементов второстепенных. Определение значимости отдельных компонентов и реакций в совокупном химизме осуществляется вычислением логарифмических откликов динамических переменных в на единичные возмущения констант скоростей. Различия в технике возмущения системных констант и способах анализа данных позволяют выделить среди этих методов основанные: 1) на функциях Грина, 2) Фурье-преобразовании, 3) стохастическом анализе, а также ряд менее распространенных вариантов. Методы функций Грина подразумевают линеаризацию уравнений, т.е. являются локальными. Применимость техники Фурье-преобразований ограничена динамической размерностью, нелинейностью и жесткостью анализируемой системы. Целью стохастического анализа ставится установить «доверительные коридоры» (см. рис. 6): их знание позволяет выявлять ведущие химические реагенты и элементарные стадии, оценивать пределы предсказуемости системной динамики на основании усеченной кинетической схемы. Следует подчеркнуть, что: 1) методы стохастического анализа во многом переплетаются с методами нелиней-
ной динамики (НД) [14]; 2) в них проявляется единство структурной устойчивости модели и ее динамики (содержания и формы; 3) получаемые с их помощью минимальные модели реакций ценны благодаря известной долгосрочности промежуточной асимптотики.
В усовершенствованной АС КИНКАТ отбор ведущих реагентов и реакционного механизма осуществляется следующим образом. Конструируется целевая функция (функционал) требуемого вида, а затем в рамках специально поставленного ВЭ отыскивается минимальное подмножество множества изначально отобранных реакций, удовлетворяющее целевому критерию. Существенно, что связи в системе (реакции) рассматриваются первичными, а элементам (реагентам) отводится второстепенная роль узлов. Опыт моделирования превращений свидетельствует (рис. 7), что в задачах ХК с N веществами число важных реакций оказывается пропорциональным ^/2. Этот факт объясняется диффузным характером межкомпонентной связности в реагирующих системах: каждое вещество способно вступать в реакции не со всеми компонентами, а только с относительно небольшим числом химически комплементарных.
Рис. 7. Номенклатурные показатели практических моделей реагирующих газов. Символом □ обозначены аналитические модели и кинетические схемы из ведущих реакций в компьютерных моделях, х - полные схемы реакций в детальных компьютерных моделях
Следовательно, в каждом элементарном акте перенос химических элементов сопровождается «локализованным» изменением фазового состояния ДС химических превращений. Ограниченность радиуса миграции компонента в фазовом пространстве за одну реакцию плюс обратимость последней означают «диффузность» переноса вещества. При данном механизме миграции удельный размер ячейки химической комплементарно-сти или числа возможных реакций, приходящихся на каждый из N компонентов, пропорционален N 1/2 Сказанное означает, что процедуру целевой редукции кинетических схем целесообразно производить в два этапа: на первом (грубом) исключаются второстепенные компоненты, а на втором (точном) - малозначительные связи между важными веществами.
Иногда кинетическую схему необходимо расширить для установления структурной устойчивости, если она хорошо описывает превращения, либо для придания динамике системы желаемых свойств в противном случае. Схема на рис. 8 иллюстрирует технику номенклатурного расширения схемы реакций применительно к горению углеводородов.
Рис. 8. Схема номенклатурного развития модели реакций в условиях подобия реагентов и элементарных актов. Свойство подобия позволяет формировать модель послойно из так называемых гомологических реакционных слоев
Как следует из представленных на рис. 8 данных, реакционные блоки второго и последующих слоев подобны базовому (слой № 1), что обеспечивает легкость их генерации на основе 1) кинетических данных для базового слоя, 2) термохимических характеристик веществ вновь формируемого слоя и 3) корреляционных методов определения кинетических констант в рядах подобных реакций. Формальным, следовательно, легко автоматизируемым способом удается составить ~ 50 % «усовершенствованной» схемы превращений: число значимых гомологических слоев п обычно составляет 2 -г- 4. Оставшиеся реакции отвечают (диффузионным в пространстве реагентов) процессам межслойного обмена: в этих процессах вещества мигрируют из слоя в слой в результате передачи друг другу гомологизующего фрагмента (например, СН2-группы). Формализовать процедуру генерации межслойных реакций сложнее, так как они не имеют прототипов в пределах базового слоя.
Отбор наиболее существенных реакций обычно основывают на физических соображениях: в удачной методике [6] таковым выступает предположение, что влияние элементарной реакции на брутто-химизм пропорционально «протекающему» через нее количеству вещества. Однако в [6] реакции рассматривались в паре с обратными: ранжирование производилось по количеству вещества, проходящего через однонаправленный канал. Очевидно, при этом «важными» оказывались все быстрые реакции, включая равновесные. С другой стороны, на динамику превращений влияет то, какие реакции быстрее -равновесные или неравновесные: ведь именно наиболее быстрые неравновесные реакции формируют химическое ядро модели. Следовательно, целевая функция (ЦФ) должна позволять «отсеять» быстрые и одновременно неуравновешенные процессы.
Эти соображения использованы авторами при формировании комплекта ЦФ в системе КИНКАТ. Пока кинетические схемы ограничивались десятком-другим реагентов и сотней реакций, «единицей отсева» при выявлении ядра химизма выбиралась реакция. При анализе современных схем превращений, в которых учтены десятки тысяч реакций между тысячами веществ и состояний, целесообразен двухэтапный подход: на грубой стадии исключаются вещества, затем на уточняющей - реакции.
Опуская несущественные детали, рассмотрим технику выделения ведущих реакций (своего рода факторного анализа) на примере чисто химических реакций. При этом отметим только, что разработанная техника применима без ограничений также в исследованиях механизма фото- и квантовохимических процессов. На грубой стадии анализа удобно контролировать динамику относительных изменений содержания отдельных химических компонентов в системе. Это осуществляется вычислением значений ff(Cj,T)lci и fb(Cj,T)lci для каждого вещества в течение всего процесса. По результатам строится динамика частотного спектра компонента. Наиболее важным компонентам отвечают широкие спектры с большими неуравновешенными амплитудами fi(Cj,T)lci Ф £ (с/,7)1с-. «Зеркальность» спектральных фрагментов ff(Cj,T)lci и fD(Cj,T)lci отвечает стационарным или равновесным состояниям. Вещества, характеризующиеся наибольшими по модулю отрицательными значениями интегрального ляпуновского показателя
Л = [ад/О]-1 (Cj,T)lCi - фШ , (20)
можно исключить из числа реагентов без ущерба для качества имитации. Переменная интегрирования выбрана в (20) таким образом, чтобы уравновесить значимость различных временных масштабов и отвечающих им групп реакций, а пределы интегрирования совпадают с длительностью эволюционного процесса в компьютерной имитации 1п/1 -г- 1п/2. На первом шаге можно отбросить половину наименее активных веществ, для которых функционал (20) дает самые дис-сипативные свойства. Затем вычисления повторяются
знсперимежалонае > ¿/теоретически — dawHöfe
Выбор реагентов Анализ
и продуктов возможных
Составление кинетической схемы
тррекиця МММ
Термокинети че ское согласование параметров ИММ, генерация программного кода
ко ррекцуя МММ
Анализ структурной устойчивости ХК-модели
ъ
Оценка динамической устойчивости ХК-модели
Рис. 9. Типовая схема составления кинетического модуля ИММ при проведении вычислительного эксперимента в газодинамических задачах технической экологии
с теми же начальными условиями, но при сокращенной вдвое номенклатуре реагентов и соответственно уменьшенном (примерно в 23/2 раза) числе реакций. Если результат усечения номенклатуры реагентов и схемы реакций (упрощенная модель химизма) удовлетворяет выбранной ЦФ, например,
ЦФ = [lnfeM)]-1 J (cN (t) - cN/2(t))/cN(t)| dint < 33 %,
i = 1, ... N/2 , (21)
процедуру повторяют. В формуле (21) cN(t) и cN/2 -соответственно кинетические кривые i-го вещества при расчете по полной и вдвое сокращенной схеме. Когда на некотором шаге оказывается, что усечение недопустимо в рамках выбранной ЦФ, в кинетическую схему «возвращается» половина наиболее значимых среди отброшенных последними компонентов. Реализующая известный метод «деления пополам» процедура выбора ведущих веществ позволяет сократить число реагентов на порядок, выполнив четыре-пять (строго говоря, log210 + 1) сравнительных расчетов на ЭВМ.
После отсева маловажных веществ можно удалить из схемы наименее важные реакции: хотя численного счета этот шаг не облегчает, он позволяет «прополоть» систему связей, избавить модель от излишнего зашумления. ЦФ возможно сконструировать самыми разными способами, выбрав в соответствии с целями исследования глобальный, как в (21), локальный или глобально-локальный критерий приемлемости. Следящими функционалами (СФ) (для простоты ограничимся здесь рассмотрением только реализованных в АС КИНКАТ СФ интегрального характера) удобно выбрать
СФ^ = [ln(t2/t1)]-1 jWf dlnt , СФиЬ = [ln(t2/t1)]-1 jWk dlnt ,
СФ2А = [lnfeM)]-1 J СФ3, = [ln(t2/ti)]-1 J
Wk - Wk Wk - Wk
dlnt ,
Wk + Wk
(22) (23)
dlnt . (24)
же на сокращение, согласно (21)-(24), при отбрасывании исключительно пар взаимно обратных реакций.
Постановка кинетической части ВЭ в технической экологии с применением АС КИНКАТ проиллюстрирована на рис. 9, алгоритма идентификации механизма превращений вынесен на рис. 10.
Функционалы (22)-(24) служат мерами: СФ^ и СФиь - относительной важности элементарных стадий, СФ2к - сравнительной важности реакционной связи в структуре минимальной модели химизма, СФ3к - степени неравновесности £-й реакции. Кандидатами на исключение из схемы являются реакции с минимальными показателями (22)-(24). Алгоритм исключения идентичен описанному выше алгоритму материальной редукции.
Подчеркнем, что, отсеивая второстепенные реакции на втором этапе сокращения кинетической схемы, мы получаем в результате существенно неравновесную минимальную модель химизма. Крайне ценное при идентификации реагирующих систем это качество делает модели, оптимизированные путем отсева элементарных процессов, принципиально непригодными в качестве элементов при синтезе больших реагирующих систем с заранее неизвестными параметрами. В последнем случае следует опираться только на результат грубой покомпонентной редукции или
Рис. 10. Схема выявления ведущих веществ и реакций в АС КИНКАТ
Работа выполнена при поддержке РФФИ (проект 05-08-33433а).
Литература
1. Моисеев Н.Н., Александров В.В., Тарко А.М. Человек и биосфера. Опыт системного анализа и эксперименты с моделями. М., 1985.
2. Варакин В.П. и др. Научно-исследовательская информационная система автоматизированного обеспечения физико-химической газодинамики. Описание проекта. МГУ, 1985. Деп. в ВИНИТИ № 2510-85.
3. Amsden A.A., O'Rourke P.J., Butler T.D. KIVA II: A computer program for chemically reactive flows with sprays / LA-11560-MS. LANL. LA, 1989.
4. Ackerman J., Wulkow M. MACRON - А program package for macromolecular kinetics // Preprint SC-90-14. Berlin, 1990.
5. Левицкий А.А., Лосев С.А., Макаров В.Н. Задачи химической кинетики в АСНИ АВОГАДРО: Мат. методы в хим. кинетике. Новосибирск, 1990. С. 7.
6. ГубановЕ.В. / Препринт / ИВТАН. № 1-344. М., 1992.
7. Давлетшин Р.Ф., Яценко О.В. // Межвед. сб.: Материалы модел. проц. управления и обработки информации. М., 1993. С. 113.
8. Губанов Е.В. и др. / Препринт / МФТИ. № 1. М., 1993.
9. Химия горения / Под ред. У. Гардинера, мл. М., 1988.
10. Tsang W., Hampson R.F. // J. Phys. Chem. Ref. Data. 1986. Vol. 15. № 3. P. 1087.
11. Барон Н.М. и др. Краткий справочник физико-химических величин. Л., 1967.
12. Основы практической теории горения / Под ред. В.В. Померанцева. Л., 1973.
13. Gear C.W. Numerical Initial Value Problems in Ordinary Differential Equations / Prentice-Hall Englewood Cliffs. N.-J., 1969.
14. Малинецкий Г.Г., Потапов А.Б. Современные проблемы нелинейной динамики. М., 2000.
Донской государственный технический университет
9 сентября 2005 г.
/
