О ВЫБОРЕ ЯЗЫКА ИНТЕРФЕЙСА ДЛЯ РЕШЕНИЯ ЗАДАЧ

4237
УДК 621.365.22-52
О ВЫБОРЕ ЯЗЫКА ИНТЕРФЕЙСА ДЛЯ
РЕШЕНИЯ ЗАДАЧ МАТЕМАТИЧЕСКОГО
МОДЕЛИРОВАНИЯ ХИМИКОТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ
С.О. Травин
ОАО ВНИИХТ
Россия, 115409, Москва, Каширское шоссе, 33
E-mail: [email protected] А.А. Быков
ОАО ВНИИХТ
Россия, 115409, Москва, Каширское шоссе, 33
E-mail: [email protected]
О.Б. Громов
ОАО ВНИИХТ
Россия, 115409, Москва, Каширское шоссе, 33
E-mail: [email protected]
П.И. Михеев
МГТУ им. Баумана
Россия, 105005, Москва, 2-я Бауманская ул., 5, стр. 1
E-mail: [email protected]
Ключевые слова: имитационное моделирование, системы ОДУ, интерфейс ввода данных, якобиан, химическая кинетика
Аннотация: Подавляющее большинство известных программных комплексов, применяемых для имитационного моделирования, используют ту или иную разновидность
языка сверхвысокого уровня, который позволяет задавать специфические для модели
вычислительные процедуры, минуя необходимость написания отдельного программного
модуля для каждой разновидности систем обыкновенных дифференциальных уравнений
(ОДУ). Не являются исключением и программные средства для решения прямой и обратной задач химической кинетики. В настоящей статье представлен подход к выбору
языка интерфейса для моделирующих комплексов химико-технологических процессов,
максимально приближенный к принятой в химии нотации для представления механизмов многостадийных реакций.
1. Введение
При рассмотрении вопросов математического и программного обеспечение для
решения прямых задач химической кинетики обычно основное внимание уделяется
подбору алгоритмов и аспектам численного интегрирования систем дифференциальных
уравнений. В то же время, нельзя не отметить, что заметное число ошибок, присутствующих в большинстве публикаций по имитационному моделированию, никоим обраXII ВСЕРОССИЙСКОЕ СОВЕЩАНИЕ ПО ПРОБЛЕМАМ УПРАВЛЕНИЯ
ВСПУ-2014
Москва 16-19 июня 2014 г.
4238
зом не связаны со степенью совершенства используемых интеграторов и сольверов, и
относятся полностью к тому, что можно назвать «человеческим фактором».
Фактически речь идет о несовершенстве программного интерфейса ввода данных
исходной задачи, который не только требует от пользователя знания более или менее
специализированного системного языка сверхвысокого уровня, но и допускает введение неточных, искаженных, ошибочных или даже внутренне противоречивых уравнений, не соответствующих химическому механизму процессов, протекающих в моделируемой системе. Казалось бы, простая задача ввода правых частей уравнений химической кинетики в случае недостаточно отработанного интерфейса ввода превращается в
труднопреодолимое препятствие. То же самое можно сказать и про задание начальных
(граничных) условий и, тем более, про ввод ограничений, задаваемых неравенствами
(например, если известна оценка сверху или снизу для диапазона значений константы
скорости какой-либо из стадий).
2. Описание алгоритма работы интерфейса ввода
Формализм прямой задачи химической кинетики достаточно прост, чтобы сделать
интерфейс ввода схемы механизма химической реакции и стартовых значений предельно простым, ясным и дружественным для пользователя, что в дальнейшем исключает возможность ошибок и/или неправильного толкования составленной системы
обыкновенных дифференциальных уравнений (ОДУ).
Непрерывное поведение химической системы на временном интервале , характеризуется вектором состояния (концентраций) ∈
и задаётся отображением
:
→
, определяющим вектор–функцию развития во времени начального состояния, при начальных условиях
.
Химическая система, содержащая веществ , ∈ 1, (включая промежуточные
нестабильные) в общем случае описывается набором химических уравнений (стадий):
→
, ∈ 1,
.
Для обратимых (или квазиравновесных) реакций отдельно учитываются прямая и
обратная стадии. Значения стехиометрических коэффициентов
участия i-го вещества в j-ой реакции (и аналогично
для продуктов реакции) a`priori считаются равными
молекулярности стадии по данному компоненту и, тем самым являются натуральными
числами в диапазоне от 1 до 3. В случае, если компонент не участвует в какой-то стадии его стехиометрические коэффициенты обращаются в ноль. В любом случае, предельное количество реагентов как в левой, так и в правой части химического уравнения
не должно превышать трех, что соответствует реальному существованию элементарных актов не выше третьего порядка:
∀
⟹ 1
3;0
3.
Различие в нижних пределах сумм стехиометрических коэффициентов для исходных реагентов (1) и продуктов реакции (0) связано с тем, что нет смысла рассматривать
процессы «самозарождения» материи из отсутствующих компонентов, тогда как превращение веществ в ходе распада в «ничто» просто соответствует отсутствию интереса
исследователя к судьбе конечных продуктов.
XII ВСЕРОССИЙСКОЕ СОВЕЩАНИЕ ПО ПРОБЛЕМАМ УПРАВЛЕНИЯ
ВСПУ-2014
Москва 16-19 июня 2014 г.
4239
В подавляющем большинстве имеющихся пакетов моделирования отображение ограничено классом систем (ОДУ) в форме Коши
, с правой частью, даваемой
законом действующих масс (для j-ой стадии):
, ∈ 1,
,
или для m-го компонента:
∙
∙
.
В компактной форме система ОДУ может быть представлена в виде:
Γ ;
где Г – результирующая стехиометрическая матрица
.
Обычно при составлении программ имитационного моделирования кинетических
систем подчинение правой части СОДУ закону действующих масс игнорируется, и
вместо этого используется гораздо более мягкий критерий Липшица. В то же время,
использование закона действующих масс не только облегчает ввод данных исключающий ошибки, но и дает возможность вычисления якобианов в явном виде, что в разы и
на порядки сокращает время работы интеграторов, использующих неявные схемы (т.е.
неразрешенные относительно ):
∙
∙
.
Тем самым выбор полиномов степени не выше третьей, соответствующих закону
действия масс, в качестве базисного набора и основного продукта для входного интерфейса оправдан не только с точки зрения надежности и удобства пользователя, но и позволяет эффективно улучшить работу интеграторов.
При таком подходе алгоритм работы интерфейса ввода выглядит следующим образом:
 Задание списка идентификаторов всех веществ, включая промежуточные, в удобной
нотации. Это могут быть привычные химические формулы (типа H2SO4 или H2O2)
или удобные мнемонические названия на понятном пользователю языке (типа серная кислота или Hydrogen_peroxide) или осмысленные идентификаторы произвольного вида. Единственным ограничением является то, что идентификатор не может
начинаться с цифры. Декларация всех значимых веществ на самой ранней стадии
позволяет избежать случайных опечаток и гарантирует от появления «двойников».
 Задание начальных концентраций. В подавляющем большинстве случаев ненулевые
значения имеют концентрации двух-трех исходных веществ. Для продуктов реакции и активных интермедиатов по умолчанию подставляются нулевые концентрации, что, впрочем, может быть легко отредактировано пользователем.
 Составляется полная схема (механизм) взаимодействия реагентов в стандартной для
химической
кинетики
нотации
(например,
→
дляреактиваФентона; здесь в составе продуктов реакции «опущен» ион гидроксония, как не влияющий на дальнейшее взаимодействие реагентов, однако если пользователя интересуют более тонкие эффекты, связанные с изменением pHсреды, он
имеет возможность включить его в состав продуктов реакции). Важно, что выбор
участников в схеме производится строго по списку, задекларированному в п. 1.
Также естественным ограничением является то, что ни среди реагентов, ни среди
продуктов элементарного акта общее количество участвующих частиц не может
превышать трех, что вытекает из требований формальной кинетики. Любые более
XII ВСЕРОССИЙСКОЕ СОВЕЩАНИЕ ПО ПРОБЛЕМАМ УПРАВЛЕНИЯ
ВСПУ-2014
Москва 16-19 июня 2014 г.
4240
сложные многочастичные взаимодействия должны конструироваться с помощью
создания промежуточных комплексов, что не составляет труда, но зато гарантирует
от использования физически бессмысленных числовых параметров.
 По умолчанию все введенные реакции имеют нулевые константы скорости. До начала расчета необходимо подставить справочные или предполагаемые значения.
Если какие-либо из стадий оставлены с нулевыми константами, их присутствие не
повлияет на результаты расчета (в том числе и на затраченное интегратором время).
Удалить реакции с нулевыми скоростями никогда не поздно, однако в ряде случаев
имеет смысл зарезервировать их с целью последующего усложнения и уточнения
модели.
Одной из первых реально работающих программных систем, основанных на перечисленных принципах, была созданная в Институте Химической Физики РАН [1]. Может показаться, что добровольный отказ от систем с более сложным видом правых частей, в том числе явно содержащих зависимость от времени, снижает функциональные
возможности системы, однако это не так. Действительно, в ряде случаев для имитационного моделирования гибридных систем [2] (например содержащих химическую часть
и управляющие воздействия программируемых клапанов) требуется задание явных
функций от времени.
В монографии [3] применительно к языку LISMA (язык инструментальных
средств) говорится: «Программное обеспечение нового приложения – задач химической кинетики со своими особенностями входного формата данных – потребовало минимальных затрат исходной спецификации входного языка».
Покажем, что инструментальный комплекс, ориентированный на «чисто химический» интерфейс, с помощью набора несложных приемов, основанных на использовании быстрореагирующего псевдовещества (при необходимости – нескольких), может
столь же просто моделировать любые дифференциальные уравнения с произвольной
правой частью. Фактически речь идет о создании вспомогательного инструментария
генерации элементарных функций с помощью чисто кинетических приемов. Арсенал
химической кинетики позволяет сделать это.
«В настоящее время имеются такого типа программы, но они, как правило, ориентированы на конкретные численные методы. Некоторые из указанных программ труднодоступны и в ряде случаев требуют дополнительного программирования со стороны
предметного пользователя. Этим обусловлена актуальность задач разработки математического и программного обеспечения для решения задач химической кинетики методом компьютерного моделирования» [4].
Самый первый шаг – это создание псевдовещества, концентрация которого будет
тождественна времени. Такое псевдовещество конструируется с помощью реакции
,
0
1.
Обратим внимание, что «генератор времени» присутствует и в левой, и в правой
частях уравнения, т.е. не расходуется. В вариантах дальнейшего использования времени в других кинетических уравнениях следует придерживаться этого же правила.
Все дальнейшие манипуляции могут проводиться как с псевдовременемtime, так и
с любыми другими веществами или псевдовеществами.
Произвольная степенная зависимость от
может быть получена
или
с помощью цепочки равновесий:
,
,
гибель,
……………………………………………………
XII ВСЕРОССИЙСКОЕ СОВЕЩАНИЕ ПО ПРОБЛЕМАМ УПРАВЛЕНИЯ
ВСПУ-2014
Москва 16-19 июня 2014 г.
4241
,
,
гибель.
Псевдовещества будут иметь концентрации равные соответствующей степени для
всех времен, превышающих τ.
Дробные степени вида 1/2, 1/3, 1/4 и т.д. могут быть получены аналогично. Например, цепочка квазиравновесий (не имеющих физического смысла) даст для вещества
p зависимость вида:
√ ,
,
гибель,
,
,
гибель,
.
Можно показать, что с помощью обычных реакций нулевого, первого и второго
порядка легко генерируются такие функции как:
 линейная функция;
 экспонента;
 гипербола.
Несколько особняком стоят такие функции как логарифм и синус, однако и они могут быть воспроизведены с помощью таких систем как брюсселятор и каталитическое
превращение с квадратичной гибелью катализатора.
В заключение приведем лишь один пример, как с помощью кинетических схем может быть построена модель истечения жидкости под собственным весом.
На рисунке представлены исходное и конечное состояние гидродинамической модели перетока жидкостей. Предполагается, что отверстие находится на высоте 20% от
верхнего края большего сосуда, соотношение площадей сечения 2:1, скорость истечения дается законом Бернулли:
1/2
Запись модели приведена ниже:
XII ВСЕРОССИЙСКОЕ СОВЕЩАНИЕ ПО ПРОБЛЕМАМ УПРАВЛЕНИЯ
ВСПУ-2014
Москва 16-19 июня 2014 г.
4242
<elementlist>
<element order="1" name="h1" description="" startvalue="100.00"/>
<element order="2" name="h2" description="" startvalue="0.00"/>
<element order="3" name="h0" description="" startvalue="80.00"/>
<element order="4" name="napor" description="" startvalue="0.00"/>
<element order="5" name="p1" description="" startvalue="100.00"/>
<element order="6" name="p2" description="" startvalue="0.00"/>
<element order="7" name="p3" description="" startvalue="0.00"/>
<element order="8" name="p4" description="" startvalue="0.00"/>
</elementlist>
<reactionlist>
<reaction order="1" source="p1+h0" product="p2" speed="10000.00"/>
<reaction order="2" source="p2" product="p1+h0" speed="1.00"/>
<reaction order="3" source="p1" product="napor+napor" speed="1.00"/>
<reaction
order="4"
source="napor+napor"
product="p1"
speed="10000.00"/>
<reaction
order="5"
source="p1+napor"
product="napor+p3"
speed="10.00"/>
<reaction order="6" source="p3+h1" product="p4+p4" speed="100.00"/>
<reaction order="7" source="p4" product="h2" speed="100.00"/>
</reactionlist>
3. Заключение
По убеждению авторов, в наилучшей степени задачам химической кинетики соответствует язык интерфейса, максимально приближенный к принятой в химии нотации
для представления механизмов многостадийных реакций. Конструирование правых
частей СОДУ, исходя из закона действующих масс, не только обеспечивает их соответствие критерию Липшица, но и дает возможность вычисления якобианов в явном виде,
что в разы и на порядки сокращает время работы интеграторов, использующих неявные
схемы.
Предлагаемый алгоритм работы интерфейса ввода содержит основные блоки:
 Задание списка идентификаторов всех веществ, включая промежуточные, в удобной нотации. Декларация всех значимых веществ на самой ранней стадии позволяет избежать случайных опечаток и гарантирует от появления «двойников».
 Задание начальных концентраций. По умолчанию для конечных продуктов реакции
и активных интермедиатов подставляются нулевые концентрации. В подавляющем
большинстве случаев «в ручную» задаются ненулевые концентрации для двух-трех
исходных веществ.
 Составляется полная схема (механизм) взаимодействия реагентов в стандартной
для химической кинетики нотации. Выбор участников в схеме производится строго
по списку, задекларированному в п. 1. Также естественным ограничением является
то, что ни среди реагентов, ни среди продуктов элементарного акта общее количество участвующих частиц не может превышать трех, что вытекает из требований
формальной кинетики.
 По умолчанию все введенные реакции имеют нулевые константы скорости. До начала расчета необходимо подставить справочные или предполагаемые значения.
Удалить реакции с нулевыми скоростями никогда не поздно, однако в ряде случаев
имеет смысл зарезервировать их с целью последующего усложнения и уточнения
модели.
Приводятся примеры решения «нехимических» моделей с использованием формального химико-кинетического языка.
XII ВСЕРОССИЙСКОЕ СОВЕЩАНИЕ ПО ПРОБЛЕМАМ УПРАВЛЕНИЯ
ВСПУ-2014
Москва 16-19 июня 2014 г.
4243
Список литературы
1.
2.
3.
4.
Брин Э.Ф., Травин С. О. Моделирование механизмов химических реакций // Хим. физика. 1991. Т.
10,.№ 6. С. 830-837.
Шорников Ю. В. Инструментально-ориентированный анализ жестких динамических, гибридных и
распределенных систем явными методами / Ю. В. Шорников, А. Ж. Абденов // ИММОД-2007: материалы 4-ой всероссийской научной конференции по имитационному моделированию. СПб., 2007.
С. 352-357.
Новиков Е.А., Шорников Ю.В. Компьютерное моделирование жестких гибридных систем. Новосибирск, 2012, 381 с.
Новиков Е.А., Шорников Ю.В. Компьютерное моделирование жестких гибридных систем. Новосибирск, 2012, 368 с.
XII ВСЕРОССИЙСКОЕ СОВЕЩАНИЕ ПО ПРОБЛЕМАМ УПРАВЛЕНИЯ
ВСПУ-2014
Москва 16-19 июня 2014 г.