CC BY
74
Вычислительные технологии
Том 11, № 1, 2006
ПРИМЕНЕНИЕ ПАКЕТА ChemPAK ПРИ МОДЕЛИРОВАНИИ ГАЗОДИНАМИЧЕСКОГО РЕАКТОРА*
В. А. Вшивков, И. Г. Черных Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: chernih@mail.ru О. П. Скляр, В.Н. Снытников Институт катализа им. Г. К. Борескова СО РАН, Новосибирск, Россия e-mail: snyt@catalysis.nsk.su
The ChemPAK program package for solving the direct problems of chemical kinetics with arbitrary number of chemical reactions is suggested. The package contains an expandable library of computational modules and provides a function of data transfer to multiprocessor. The process of Ci — C2 hydrocarbon's pyrolysis in gas-dynamics reactor with emission has been numerically simulated using the ChemPak software package.
Введение
Математическое моделирование в области физико-химической газодинамики реагирующих сред широко используется при создании новых и усовершенствовании имеющихся способов переработки добываемого углеводородного сырья. Для переработки природного газа, состоящего в основном из углеводородов, перспективны процессы пиролиза в газодинамических реакторах с получением более тяжелых соединений [1]. При математическом моделировании процессов пиролиза необходимо решать систему уравнений динамики газа вместе с уравнениями химической кинетики. Последние уравнения представляют собой системы обыкновенных дифференциальных уравнений (ОДУ), описывающих кинетические превращения отдельных компонентов реакционного газа, число которых достигает сотен и более с соответствующей размерностью системы ОДУ. Для решения задач химической кинетики с большой размерностью необходимо современное программное обеспечение, удовлетворяющее ряду требований. Среди них следует указать наличие эргономичного интерфейса и расширяемой библиотеки вычислительных методов, возможность работы с современными базами физико-химических данных (GRIMECH [2], NIST [3], NASA [4] и др.), возможность работы в связке персональный компьютер — суперЭВМ. В последнее
*Работа выполнена при финансовой поддержке интеграционных проектов СО РАН № 122 и № 148, программы Президиума РАН № 25-2, Российского фонда фундаментальных исследований (грант № 0501-00665).
( Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
время появился ряд программных продуктов, ориентированных на описанные выше задачи. К крупным программным продуктам для широкого круга задач моделирования можно отнести FLUENT [5], CHEMKIN [6], StarCD [7], HYSYS [8] и др. Существуют и небольшие программные пакеты (CKS [9], Kintecus [10], AcuChem [11], ChemMathS [12] и др.). Кроме указанных коммерчески распространяемых пакетов программ существует также ряд специализированных библиотек подпрограмм NAG [13], Numerical recipes [14], БанКин [15], разработанных в разные годы, и компьютеризированный справочник "Физико-химические процессы в газовой динамике" [16]. К недостаткам крупных пакетов можно отнести высокую цену и, как правило, отсутствие возможности расширения и усовершенствования пакета пользователем. Недостатками специализированных библиотек подпрограмм являются высокая сложность структуры данных, отсутствие эргономичного интерфейса, а также сложности, связанные с архитектурой библиотек, не рассчитанной на использование ее в комплекте с единым интерфейсом. Чтобы преодолеть указанные недостатки, в работе [17] предложена библиотека классов для решения прямых задач химической кинетики. В ней продемонстрирован объектно-ориентированный подход к созданию пакета для расчета термодинамических свойств и химической кинетики для химических реакций. Однако использование этой библиотеки для решения задач химической кинетики вместе с другими уравнениями требует создания дополнительных программных интерфейсов, в том числе для работы в связке с суперЭВМ.
Таким образом, проведенный анализ имеющихся программных пакетов и библиотек программ показал необходимость разработки нового программного пакета, ориентированного на решение прямых задач химической кинетики в сетевой среде с параллельными суперЭВМ. В нем должны быть решены следующие задачи:
— ввод и изменение системы уравнений большого (свыше 1000) числа химических реакций в общепринятой в химии нотации;
— проверка корректности записанной системы уравнений химических реакций;
— нахождение для заданной системы уравнений химических реакций правых частей ОДУ, которые должны быть записаны в универсальной нотации алгоритмов, работающих на персональных компьютерах и параллельных суперЭВМ;
— создание итеративного режима работы с возможностью проверки решений укороченных систем ОДУ на однопроцессорных ЭВМ;
— представление полученного блока решения системы ОДУ для химических реакций в виде программного модуля, встраиваемого в программы решения систем уравнений математической физики, моделирующих на ЭВМ изучаемые процессы.
Созданный программный пакет ChemPAK решает поставленную задачу. Ниже приведено описание принципов его работы и продемонстрировано использование на задаче математического моделирования процессов в газодинамическом реакторе с излучением.
1. Программный пакет ChemPAK
Прямые задачи химической кинетики решаются исследователем посредством задания химических реакций для реагентов, получения систем ОДУ для реакций, решения этих систем. Аналогично для пользователя решение задач с пакетом ChemPAK было разделено на три основных этапа. На первом из них происходит работа с базой данных — ввод и редактирование химических реакций, а также различных химических данных, необходимых для моделирования. На втором — работа с транслятором (трансляция системы химических
Рис. 1. Схема программного пакета.
реакций в систему обыкновенных дифференциальных уравнений, генерирование матрицы Якоби для полученной системы, расчет скоростей протекания химических реакций и др.). Третий этап — работа с вычислительными модулями, в частности решение системы уравнений, описывающей химическую кинетику, совместно с уравнениями газодинамики. Исходя из такого принципа работы созданный программный пакет был функционально разделен на три части (рис. 1): база данных химических реакций, химический транслятор, вычислительные модули.
1.1. База данных химических реакций
База данных химических реакций создана для хранения систем химических реакций, а также вспомогательной информации, необходимой для генерации систем ОДУ и численного моделирования. Для работы с базой данных разработан отдельный модуль, являющийся частью транслятора химических реакций.
Для разработки базы данных использовалась СУБД Interbase. Наличие бесплатной, свободно распространяемой версии, скромные требования к ресурсам компьютера (для работы достаточно процессора Pentium 2 и 64 Мбайт оперативной памяти), а также наличие библиотек прямого доступа (без использования интерфейсов типа ODBC, ADO) к базе данных обусловили выбор СУБД. К тому же Interbase имеет высокую производительность и удобный интерфейс для работы с приложениями, написанными в Borland C+—+ Builder. Для обеспечения возможности работы нескольких пользователей с одним набором данных доступ к базе данных реализован через локальную или глобальную сеть.
Формат хранения данных в виде разобранной по химическим элементам реакции был выбран из соображения увеличения скорости трансляции системы химических реакций. Хранение данных в таком виде несколько увеличило размер базы данных. Но этот подход позволил значительно увеличить скорость трансляции благодаря отсутствию необходимости делать разбор строки с химической реакцией перед транслированием. Для увеличения
скорости трансляции также разработан алгоритм кэширования данных, который позволяет минимизировать обращения к базе данных в процессе трансляции.
На рис. 2 приведена логическая модель представления химических данных в базе (стандарт IDEF1X). Модель описывает химические данные как факты о сущностях и связях между сущностями. Кратко методология IDEF1X изложена в работе [18]. На рис. 2 прямоугольниками обозначены сущности, в частности алфавит химических элементов, химические реакции, вспомогательные данные. Сплошными и пунктирными линиями показаны связи между сущностями с краткими текстовыми комментариями, поясняющими тип связи. Идентифицирующая связь между сущностью-родителем и сущностью-потомком изображается сплошной линией. Сущность-родитель в идентифицирующей связи может быть как независимой, так и зависимой от идентификатора (это определяется ее связями с другими сущностями). Пунктирная линия изображает неидентифицирующую связь. Сущность-потомок в неидентифицирующей связи будет независимой от идентификатора, если она не является также сущностью-потомком в какой-либо идентифицирующей связи. Моделирование базы данных производилось с помощью пакета Allfusion ERWin Data Modeller [19].
Справочник таблиц химических реакций
Алфавит химических символов
Alphabet
Реакции
Вспомогательные данные состоят из
Table Name Backup Exist Heat Table Exist
Substance Capacity
Константы используют
Наличие таблиц теплот
Состоит из
Таблица химических реакций
Reaction Number Reaction Reaction Type
Production Velocity Coeff Direct Production Velocity Coeff Reverse Src Literature
Наличие резервной
Таблица теплог
-'-4-
Reaction Number Chemical Element
Резервная таблица химических реакций
Reaction Number Reaction Reaction Type
Production Velocity Coeff Direct Production Velocity Coeff Reverse Src Literature
Рис. 2. Схема базы данных химических реакций в пакете СЬешРЛК.
Преимуществами данного пакета являются возможность генерирования готового SQL-кода базы данных, а также обратное проектирование модели. Благодаря обратному проектированию можно внести изменения в базу данных, и затем ERWin Data Modeller сам внесет изменения в модель. Используя эту возможность, удобно документируется база данных и наглядно проверяется корректность связей между сущностями, созданными в процессе работы базы данных.
1.2. Транслятор
Современные концепции создания многомодульных систем (пакетов программ) предполагают использование единого интерфейса для всех модулей программного пакета, а также возможность встраивания новых модулей без необходимости новой сборки ядра. Химический транслятор отвечает этим требованиям и является ядром программного пакета. Для оптимизации процесса разработки архитектуры транслятора было использовано CASE-средство Rational Rose. Благодаря возможности в короткий срок создавать многооконные интерфейсы любого уровня сложности, наличию библиотеки VCL, облегчающей реализацию алгоритмов обработки текстовых данных, наличию библиотеки IBX (библиотеки прямого доступа к базам данных СУБД Interbase) для создания транслятора была выбрана среда разработки Borland C+—+ Builder. Особенности представления данных в задачах химической кинетики и алгоритм трансляции систем химических реакций [20] из классического вида в систему ОДУ предполагают использование специализированных типов
Рис. 3. Функциональная модель транслятора.
данных, необходимых для эффективной работы транслятора. Для ускорения работы алгоритмов трансляции система химических реакций кэшируется из базы и хранится в оперативной памяти. Для хранения в памяти исходной и оттранслированной системы реакции используется класс AnsiString библиотеки VCL. Это дает отсутствие ограничений на длину строк (она ограничена лишь объемом оперативной памяти), эффективное использование памяти, набор методов работы со строками (поиск подстроки в строке, сравнение строк, замена части строки на другую строку, удаление части строки, склейка строк), реализованных на языке Ассемблер. На рис. 3 представлена функциональная модель транслятора.
Для подключения вычислительных модулей без необходимости новой сборки пакета был разработан интерфейс работы с динамически подключаемыми библиотеками (DLL). Формат выходных данных транслятора создан с учетом возможности работы со старыми библиотеками программ и вычислительными модулями, находящимися на многопроцессорных компьютерах. Благодаря этому решена задача создания универсальной нотации для хранения оттранслированных данных.
Ниже в качестве примера на языке Фортран приведены фрагменты выходных файлов с системой ОДУ и вспомогательными данными, соответствующие системе химических реакций этилена.
Файл etilen_pyrolysys.dat (Файл содержит систему ОДУ) f(1) =
&-k(1)*y(1)
&*y(2)
Файл etilen_pyrolysys_map.dat (Файл содержит химические элементы, участвующие в системе химических реакций) y(1)=C2H4 y(2)=M y(3)=C2H2 y(4)=H2 y(5)=C2H3 y(6)=H
Файл etilen_pyrolysys_k_T_1200.dat (Файл содержит константы скоростей протекания реакций)
k(1) = 1228,15454101563 k0(1)=2187181824
Для работы с широким кругом химических данных, представленных в различных источниках (статьях, базах данных различных разработчиков), создан модуль обмена данными с приложением Microsoft Excel. Большинство данных из баз сторонних производителей представимо в формате таблиц Microsoft Excel. C помощью аппаратно-программных средств сканирования текста можно системы химических реакций из печатных источников перевести в Excel и в дальнейшем — в единую базу данных. Для импортирования и экспортирования системы химических реакций с константами скоростей протекания используется библиотека Excel_2K_SRVR. Выбор этой библиотеки обусловлен простотой в использовании и широким набором методов работы с приложением Microsoft Excel. Ниже приведен фрагмент программы, реализующий импорт констант скоростей протекания
химических реакций из таблицы Excel в базу данных химических реакций. Аналогичным способом реализован экспорт данных. try
{
Ехсе1Арр11са^оп1->Соппес^);//Подключение транслятора к Microsoft Excel }
catch (...) {
MessageD1g("Microsoft Excel не установлен", mtError, TMsgD1gButtons() « mbYes, 0);//Сообщение об ошибке в случае невозможности подключения
Abort; }
for (int i = StrToInt(RangeStart->Text);i<=StrToInt(RangeEnd->Text);i++)
//Основной цикл считывания данных из Excel {
Exce1App1ication1->Workbooks->get_Item(ItemIndex);//получение ссылки на файл Excel
Exce1Workbook1->ConnectTo(Exce1App1ication1->Workbooks->get_Item(ItemIndex)); //Подключение к файлу Exce1
Exce1Worksheet1->ConnectTo(Exce1Workbook1->Worksheets->get_Item(ItemIndex)); //Подключение к странице файла с данными Excel RangePtr R;
R = Exce1Worksheet1->get_Range(TVariant(Exce1Co1umn->Text+IntToStr(i)), TVariant(Exce1Co1umn->Text+IntToStr(i)));//Получение ссылки на элемент столбца с данными
TVariant abc = R->get_Va1ue();//Получение значения
tempVal = F1oatToStrF(abc.db1Va1,ffExponent,15,2);//Преобразования формата числа Exce1 во внутренний формат базы данных химических реакций MainForm->IBQuery2->Active = false;
if (DirectWay->Checked) MainForm->IBQuery2->SQL->Text = "UPDATE " +MainForm->Tab1eNum->Text.UpperCase()+"SET K_NUM = '"+tempVa1+ \ WHERE REACTION_N = "+IntToStr(i);
else MainForm->IBQuery2->SQL->Text = "UPDATE "+
MainForm->Tab1eNum->Text.UpperCase()+"SET K2_NUM = '"+tempVa1+ \ WHERE REACTION_N = "+IntToStr(i);//Составление запроса на добавление импортированных данных
MainForm->IBQuery2->Open();//Выполнение запроса на добавление импортированных
данных }
MainForm->IBTransaction->CommitRetaining();//Подтверждение добавления импортированных данных Exce1App1ication1->Connect();
Как видно из фрагмента программы, работа с данными таблиц Excel реализована через три объекта — ExcelApplication, ExcelWokrbook, ExcelWorkSheet. Обращение к этим объектам реализовано через стандартный механизм get/set методов. Сам механизм стыковки приложений скрыт в библиотеке Excel_2K_SRVR. В свою очередь, благодаря мощному интерфейсу Microsoft Excel можно использовать для работы с различными химиче-
скими данными. Например, широко используемая база данных GRIMECH [2] через web-интерфейс базы данных быстро переносится в Excel обычным копированием web-страницы и через импорт переносится в базу данных ChemPAK. Excel также использовался для конвертации большого объема числовых данных из одного формата в другой. Например, конвертация 200 чисел из формата ккал/моль в формат кДж/моль с экспортом и импортом заняла несколько секунд.
1.3. Использование транслятора химических реакций
С точки зрения пользователя работа с транслятором химических реакций производится следующим образом.
Этап 1 — ввод системы химических реакций с автоматической проверкой правильности написания химических реакций или импорт системы из сторонней базы с аналогичной проверкой правильности написания химических реакций.
Этап 2 — ввод констант скоростей протекания реакций (покомпонентный ручной ввод, ручной ввод уже рассчитанных констант, импорт рассчитанных констант из Microsoft Excel).
Этап 3 — запуск трансляции системы химических реакций в систему ОДУ.
Этап 4 — вычисление матрицы Якоби для системы ОДУ (необходимо использовать только в том случае, если матрица используется в численном решении ОДУ).
Этап 5 — расчет скоростей протекания химических реакций в случае изотермичности процесса (расчет проводится только в том случае, если данные для констант заданы в развернутом виде. Если они уже рассчитаны, то происходит вывод констант в текстовый файл напрямую из базы данных).
Этап 6 — расчет и вывод вспомогательных данных (теплоты образования радикалов в различных форматах).
1.4. Вычислительные модули
Вычислительные модули программного пакета предназначены для численного решения полученных систем ОДУ для записанной схемы химических реакций при задаваемых начальных данных. Пользователь имеет возможность воспользоваться как включенными в состав пакета методами, так и своими собственными процедурами исходя из особенностей решаемой задачи. В программный пакет для жестких систем включены диагонально-неявный метод Рунге — Кутты RADAU5, имеющий пятый порядок аппроксимации, и полунеявный метод Розенброка четвертого порядка RODAS. Подробное описание особенностей реализации этих методов и результатов численных экспериментов приведено в [22]. Проведенные нами расчеты для систем химических реакций, включающих несколько десятков компонентов, подтвердили вывод [22], что использование RADAU5 целесообразно в тех случаях, когда требуемая точность расчетов выше 10-4. В случаях, когда такая точность не требуется, эффективно применяется RODAS. Для решения нежестких задач применяется процедура DOPRI8, созданная авторами [21] и включенная в пакет ChemPAK. Численный алгоритм основан на явной схеме Дорманда и Принса восьмого (седьмого) порядка аппроксимации с процедурой выбора шага по вложенной формуле. Эти методы в свободно распространяемой реализации [21, 22] в настоящее время широко используются
для интегрирования систем ОДУ. При необходимости работы с внешними вычислительными модулями их подключение осуществляется динамически в виде библиотеки DLL.
Код сгенерированных ChemPAK может быть включен в программы численного моделирования физико-химических процессов как на персональных, так и на суперЭВМ, в том числе Fluent [5].
2. Результаты тестирования транслятора
Для тестирования транслятора использовалась следующая конфигурация компьютера: AMD Athlon XP 2600 (1.9 ГГц), RAM 768 Мбайт. Тестирование производительности транслятора проводилось на системах из 200, 400, 800, 2000, 4000, 8000 модельных реакций. На рис. 4 приведены два графика зависимости времени трансляции от количества химических реакций (штриховая линия — трансляция химических реакций и генерация матрицы Якоби, сплошная линия — только трансляция химических реакций).
Видно, что эффективность алгоритма трансляции приблизительно порядка N. Количество уравнений системы ОДУ, соответствующей системе химических реакций, зависит от количества химических веществ, участвующих в реакциях. Обычно для задач пиролиза углеводородов количество химических веществ и радикалов в 5-10 раз меньше количества химических реакций. Алгоритм трансляции реализован таким образом, что одно уравнение получается за один проход всей системы химических реакций. Благодаря этому достигнута высокая эффективность алгоритма трансляции. Генерирование матрицы Якоби не оказывает существенного влияния на производительность транслятора. Кэширование данных из базы данных также позволило увеличить скорость трансляции.
35
2000
4000
Количество реакций
6000
8000
Рис. 4. Производительность транслятора.
3. Математическое моделирование газодинамического реактора с излучением
3.1. Постановка задачи
Примером использования созданного пакета СЬешРАК служит изучение динамики реагирующего потока в коническом сопле с дополнительным воздействием излучения. К целям математического моделирования относилось исследование различных предполагаемых кинетических схем разложения углеводородов, а также изучение влияния газодинамических параметров и начальных данных (расхода, мощности излучения, состава газовой смеси) на режим работы реактора.
Схема исследуемого процесса приведена на рис. 5. Газообразный реагент движется с дозвуковой скоростью в сужающейся трубе. В трубу вводится энергия в виде излучения, которая имеет большое значение для протекания химических реакций.
Для описания процесса в упрощенном случае (с зависимостью только от продольной координаты г) предложена следующая математическая модель:
а
-Т (Яр*) = 0,
аг
а ая
т(я ^+Р»=^
( ат а
РЧ+ (7 - 1)таг (Яу)
— (/Я) = —аоП24/Я, аг
—а (/я ) — 2
Як(Т)
До
1 --)
о
(Т - ) + £
— (ЯнгУ) = г Е 1,3;
а
0,
Ао = ЯоРоУо - расход,
Т = То — начальная температура,
/ = /о — плотность мощности излучения,
н
н
— объемные доли реагентов в газовой смеси,
¿1, р = р 1
давление на выходе.
Для переменных величин и параметров системы использованы следующие обозначения: Я — площадь поперечного сечения трубы; у — среднемассовая скорость; р — плотность; р — давление газа; нг — объемная концентрация г-го вещества; / — плотность
г
о
2
Поток газа
Рис. 5. Схема реактора.
мощности излучения; п24 — объемная концентрация С2Н4; — молекулярная масса ¿-го вещества; Т — температура газовой смеси; Тш — температура стенки сопла; а0 — коэффициент поглощения излучения; /шг — скорость изменения концентрации ¿-го вещества, вычисленная по закону действующих масс из схемы химических реакций; ^ — теплота образования (распада) ¿-го вещества; 7 — показатель политропы; ср — теплоемкость смеси при постоянном давлении; к' — коэффициент теплопроводности для границы горячего газа с холодным; Ко — начальный радиус трубы; — координата вершины конуса; К — универсальная газовая постоянная.
Здесь 3 — количество веществ, участвующих в химической реакции:
Б(г) = 1 -
Здесь 50 = пК0 — площадь входного сечения.
Плотность газа связана с концентрациями реагентов соотношением р ^ ^ цгпг. В ка-
г
К т
честве уравнения состояния используется р = — рТ .
ц
Три первых уравнения представляют собой уравнения газовой динамики (неразрывности, движения, энергии). Система ОДУ включает также уравнение для излучения (закон Бугера — Ламберта) и уравнения химической кинетики, описывающие изменения концентраций веществ в реагирующей смеси. Уравнение для температуры газа дополнено слагаемыми в правой части, связанными с выделением энергии от излучения, тепловыми потоками на стенки сопла и тепловыми эффектами химических реакций.
В качестве граничных условий при г = 0 задаются расход, температура газа, плотность мощности излучения и доли реагентов в составе газовой смеси. На правой границе задается давление реакционной смеси. Тем самым решается краевая задача. Давление и абсолютные концентрации химических элементов на левом сечении реактора определяются решением.
3.2. Применение СЬешРЛК
В настоящее время кинетическая схема химических реакций неразбавленных смесей углеводородов является предметом интенсивных исследований. Для изучения динамики реагирующего потока были отобраны три кинетические схемы пиролиза природного газа и этилена с разным числом реакций, включая традиционно применяемую схему Касселя [25] и более детальные схемы, предложенные в последнее время для диапазона температур смеси, в котором должен работать изучаемый реактор [26, 27]. Используя СЬешРАК, системы химических реакций и данные для расчета скоростей реакций и теплот образования (распада) реагентов ввели в базу данных. Соответствующие блоки из скоростей образования-расходования реагентов /шг сгенерировали транслятором в виде фортран-кода. Систему уравнений газовой динамики записали в переменных р, I, Т, пг и разрешили относительно производных, после чего в нее включили блоки /шг.
Для нахождения численного решения краевой задачи был реализован метод пристрелки, использующий включенные в состав пакета модули.
3.3. Результаты моделирования
Химические процессы в реакторе описываются кинетическими схемами превращения углеводородов с большим числом неизвестных параметров, переменными условиями нагрева газовой смеси излучением. На химические реакции основное влияние оказывают температура, состав и скорость смеси. Кинетическая схема должна удовлетворительно передать соотношение реагентов, фиксируемых газовым хроматографом (Н2, С1— С4), и дать прогноз на появление более тяжелых углеводородов. В число реагентов, фигурирующих в кинетических схемах, входили водород и углеводороды, относительная молекулярная масса которых колебалась от 14 (СН2) до 172 (С14Н4).
3.3.1. Схема Касселя
Одной из наиболее простых и распространенных для первичного анализа схем разложения метана является схема Касселя [25] с параметрами:
2СН4 —> 1 С2Н4 + 2Н2,
С2Н4
С2Н2
М
С2Н2 + Н2,
м
2С + Н2,
к1 = 5.0 ■ 1013 ■ ехр к2 = 9.0 ■ 1013 ■ ехр
ЕаЛ 1
КГ) с,
к3 = 1.78 ■ 104
ехр
Еа2 КТ
Еаз КТ
1
с 1
с
Еа1 = 87
Еа2 = 69 Еа3 = 23
ккал моль ккал моль ккал моль
В схеме отсутствуют С3 и более тяжелые углеводороды. Для получения численных решений была использована процедура И,АВАи5, поскольку данная кинетическая схема приводит к жесткой системе ОДУ.
На рис. 6 представлены сравнительные результаты расчета режима реактора без учета химических реакций (штриховая линия) и с механизмом протекания реакций по схеме Касселя (сплошная линия) при следующем наборе параметров и начальных данных: К0 = 6мм, ¿0 = 100 мм, Тш = 600 К, Т0 = 500 К, А0 = 12л/ч, пС2н4 = 20%, пСН4 = 80%,
7 = 5/3, Ср =
-1
р1 = 1 атм, К = 5 Дж/(К ■ м2 ■ с), 10 = 12 ■ 105 Дж/(м2 ■ с), а0 = 133.5м 2150м2/(с2 ■ К).
Графики приводятся для безразмерных переменных (температуры, концентраций метана и этилена, отнесенных к концентрации газовой смеси, и плотности мощности излучения). Температура газового потока обезразмерена на температуру стенки, плотность мощности излучения — на начальную плотность мощности.
Как следует из графиков, учет химических реакций оказывает существенное влияние на распределение газодинамических параметров системы. С уменьшением концентрации этилена в составе смеси снижается поглощение энергии излучения, что совместно с эндотермическим характером химических реакций приводит к падению температуры смеси в сравнении со случаем отсутствия реакций.
Результаты расчетов показывают, что реакционной зоне, где происходит основное изменение концентраций метана и этилена, предшествует зона нагрева газа под действием
Рис. 6. Влияние химических реакций на параметры системы: сплошная линия — с учетом химических реакций, штриховая — без него.
излучения. При указанных параметрах и начальных данных реакционная зона формируется в ограниченной области по достижении температурой газовой смеси значения около 1200 К. Сама реакционная зона имеет длину примерно 0.35... 0.4 мм. Внутри нее температура смеси меняется незначительно.
3.3.2. Схема дегидроконденсации метана в неравновесных условиях
Более подробный механизм дегидроконденсации метана с соответствующими параметрами реакций был предложен в работе [26] для экспериментов с ударными волнами. Он состоит из 42 реакций для 19 реагентов. Но нам не удалось по этой кинетической схеме рассчитать химические превращения во всем газодинамическом реакторе с указанными плотностями мощности излучения. Как оказалось, данная кинетическая схема применима в узком относительно заданных условий диапазоне параметров по давлению, температуре и составу смеси. Попытка выйти за этот диапазон приводила к неприемлемо большой жесткости системы ОДУ. Работа с программным пакетом СЬешРАК в итеративном режиме с быстрым внесением изменений в кинетическую схему и в набор термодинамических параметров приводила к созданию новых кинетических схем, отбор среди которых может быть проведен лишь с привлечением экспериментальных данных. Этот пример показал необходимость дополнить пакет СЬешРАК программой с алгоритмом сокращения числа химических реакций для снижения жесткости получаемой системы ОДУ.
3.3.3. Схема пиролиза метан-этиленовой смеси
Схема кинетических реакций пиролиза этилена с участием метана предложена в работе [27] для широкого диапазона температур. Здесь же приводятся кинетические константы в виде зависимостей, как и в схеме Касселя, для расчета скоростей протекания химических реакций. Эта схема включает 85 прямых и столько же обратных реакций для 30 хими-
1.6 -г
Температура
Температура
О 0.25 0.5
0.1 -г Концентрация С2Н4
1 \
0.05 --
0.25 0.5
0.75 1
г
0.15 у Концентрация С2Н4
0.1
0.05
0.25 0.5
0.75 1
2
\
0.9 - V
0.8 --
0.7
Плотность мощности излучения
0.25 0.5 0.75 1
2
Сплошная линия — расход 30 л/ч. Штриховая линия — расход 50 л/ч. Штрихпунктирная линия — расход 100 л/ч. Смесь содержит 10 % этилена.
0.95
Плотность мощности излучения
Сплошная линия — 5 % этилена. Штриховая линия — 10 % этилена. Штрихпунктирная линия — 15% этилена. Расход смеси 30 л/ч.
Рис. 7. Влияние расхода и состава газовой смеси на параметры системы.
ческих веществ, протекающих в газовой фазе. Расчеты реактора с такой схемой методом RADAU5 в зависимости от задаваемой точности решения занимают на Pentium 4 CPU 2.6 ГГц от 3 до 20 с. Ясно, что решение двух- и трехмерных задач газовой динамики с учетом сложных кинетических эффектов потребует применения параллельных суперЭВМ.
На рис. 7 показано влияние расхода и количества этилена в составе смеси на динамику процесса. В качестве базового взят расчет при следующих параметрах и начальных данных: R0 = 6 мм, z0 = 100 мм, Tw = 600 К, T0 = 500 К, A0 = 30 л/ч, Пс2н4 = 10%, nCH4 = 90%, p1 = 1 атм, К = 5 Дж/(К ■ м2 ■ с), 10 = 12 ■ 105 Дж/(м2 ■ с), а0 = 130.6м-1, 7 = 5/3, Cp = 2150м2/(с2 ■ К).
Из рис. 7, a и б следует, что чем выше расход газа, тем дальше от начального сечения сопла располагается реакционная зона. При изменении расхода и состава смеси длина реакционной зоны, имея значение порядка 0.1, меняется незначительно и зависит от кинетических параметров схемы. Кинетические данные этой схемы (в отличие от схемы Касселя) обеспечивают полный распад этилена в зоне реакции. С окончанием распада этилена передача энергии излучения в среду и последующий нагрев газовой смеси не осуществляются. В результате плотность мощности излучения остается постоянной. Зависимость температуры, при которой начинаются реакции по исследуемой кинетической схеме, от начального количества этилена в составе газовой смеси показана на рис. 7, г. Стартовая температура колеблется от 840 (при 5% этилена) до 900K (при 15% этилена).
3.3.4. Обсуждение результатов
На рис. 8 представлены результаты моделирования для двух кинетических схем при указанных в подразд. 3.3.3 параметрах и начальных данных. Как следует из графиков, разные схемы приводят к сильно отличающимся расчетным режимам. Причину расхождения результатов, полученных по схеме Касселя и схеме [27], следует искать в применимости этих схем для описания экспериментальных данных. Ясно, что обе эти схемы не вполне пригодны для моделирования пиролиза легких углеводородов в реакторе: схема Касселя не учитывает появления более тяжелых соединений, а схема пиролиза этилена, по-видимому, содержит кинетические параметры для реакций углеводородов, не соответствующие изучаемым условиям. Решение этих проблем лежит на пути конструирования новых кинетических схем на базе наиболее простых с применением данных по механизмам реакций
Концентрация С2Н4 N
\
\
\
Ч
.25 0.5 0.75 1
z
Температура
0.1 п
0.075 -
0.05 -
0.025 -
0 0
Рис. 8. Результаты расчетов по двум кинетическим схемам: сплошная линия — схема пиролиза этилена, штриховая — схема Касселя.
и скоростям их протекания из различных источников. Используя пакет ChemPAK, можно оперативно оценивать применимость литературных данных по отдельным стадиям и в целых кинетических схемах для изучаемых экспериментальных условий. С его помощью легко конструируются гипотетические схемы и проводятся тестовые расчеты прямых задач химической кинетики.
Установить наиболее адекватную кинетическую схему можно путем сопоставления экспериментальных данных с расчетными.
Выводы
Для решения прямых задач химической кинетики создан программный пакет ChemPAK с расширяемой библиотекой вычислительных модулей, позволяющий обрабатывать большие системы химических реакций, записанных традиционными химическими символами. Пакет ориентирован для работы на однопроцессорных ЭВМ с возможностью передачи подготовленных данных на суперЭВМ. Одним из назначений пакета является исследование, посвященное определению механизмов химических реакций, протекающих в разных газодинамических условиях.
С помощью представленного пакета исследована проблема применимости существующих представлений о механизме пиролиза легких углеводородов в газодинамическом реакторе с вводом энергии посредством излучения.
Список литературы
[1] Буянов Р.А., Васильева Н.А., Пармон В.Н. и др. Эндотермический химический реактор с газодинамическим управлением. Новосибирск, 2001 (Препр. РАН. Сиб. отд-ние. ИТПМ. № 5-2001).
[2] http://www.me.berkeley.edu/gri_mech/
[3] http://www.nist.gov/
[4] http://www.nasa.gov/
[5] http://www.fluent.com/
[6] http://www.ca.sandia.gov/chemkin/
[7] http://www.cd-adapco.com/
[8] http://www.aspentech.ru/RU/products/eng/hysys/process/
[9] https://www.almaden.ibm.com/st/computational_science/ck/msim/
[10] http://www.kintecus.com
[11] http://www.enveng.ufl.edu/homepp/andino/Extra_2002.htm
[12] http://www.cesd.com/cesddls.html
[13] http://www.nag.co.uk/
[14] http://www.nr.com/
[1Б] Бабкин B.C., Бабушок В.И., Дробышев Ю.П. и др. Автоматический банк кинетической информации, общее описание. 19В7 (Препр. ВЦ СО АН СССР. Т. 704).
[16] Компьютеризированный справочник "Физико-химические процессы в газовой динамике" / Под ред. Г.Г. Черного и С.А. Лосева. Т. 1. Динамика физико-химических процессов в газе и плазме. М.: Изд-во МГУ, 199Б.
[17] Мигун А.Н., Матвейчик Е.А., Чернухо А.П., Жданок С.А. Объектно-ориентированный подход к моделированию химической кинетики // ИФЖ. 200Б. Т. 78, № 1. С. 153-158.
[1В] http://www.citforum.ru/database/case/glava2_4_2.shtml
[19] http://www.ca.com/
[20] Замараев К.И. Химическая кинетика: Курс лекций. Новосибирск: Изд-во НГУ, 1994.
[21] ХаЙрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений, нежесткие задачи. М.: Мир, 1990.
[22] ХаЙрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений, жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999.
[23] Kuptsov P.V., Kuznetsov S.P., Knudsen C. Connective wave front for a reactor-diffusion system in a conical flow reactor // Phys. Letter. 2002. Vol. 294. P. 210-21б.
[24] Тынников Ю.Г., Генкин В.Н. Газокинетическая схема высокоэнергетического воздействия на поток метана // Письма в ЖТФ. 2000. T. 24, № 24. С. б4-б9.
[25] Артамонов А.Г., Володин В.М., Авдеев В.Г. Математическое моделирование и оптимизация плазмохимических процессов. М.: Химия, 1989.
[26] Yoshiaki Hidaka, Kazutako Sato, Yusuke Henmi et al. Shock-tube and modelling study of methane pyrolysis and oxidation // Combustion and Flame. 1999. Vol. 118. P. 340-358.
[27] ЖильцовА И.В., ЗАСлонко И.С., Карасевич Ю.Л., Вагнер Х.Г. Неизотермические эффекты в процессе сажеобразования при пиролизе этилена за ударными волнами // Кинетика и катализ. 2000. T. 41, № 1. С. 87-101.
Поступила в редакцию 31 мая 2005 г., в переработанном виде —12 октября 2005 г.
