CC BY
72
УДК 544.45
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ АДИАБАТИЧЕСКОГО ПЕРИОДА ИНДУКЦИИ ДЛЯ МЕТАН-КИСЛОРОДНЫХ СМЕСЕЙ В ШИРОКОМ ДИАПАЗОНЕ НАЧАЛЬНЫХ ДАВЛЕНИЙ И ТЕМПЕРАТУР
В.К. Рябинин, Ю.М. Ковалев
Исследуются особенности процессов моделирования теплового взрыва в изохорном адиабатическом реакторе в широком диапазоне начальных давлений и температур для метан-кислородных и метан-воздушных смесей с использованием максимальных кинетических механизмов разветвленно-цепных реакций. Рассмотрены методы построения математических моделей на базе уравнений состояния идеального и реального газов. Учтены поправки к энтальпии и теплопроводности компонентов газовой смеси за счет высокого давления в реальном газе. Сравниваются и анализируются результаты расчетов для этих двух классов моделей. Показано, что при высоких давлениях период индукции в реальном газе больше, чем в идеальном при тех же начальных условиях. Установлено, что доминирующим фактором, определяющим этот эффект, является разница в плотности идеального и реального газа после начального сжатия.
Ключевые слова: математическое моделирование, горение и взрыв, цепные реакции, высокие давления.
Оценка величины адиабатического периода индукции и (часто называемого временем задержки воспламенения) для реагирующих смесей произвольного состава при заданной начальной температуре и давлении необходима для оценки пожаро-взрывобезопаспости веществ и технологических процессов, а так же для современных методов анализа и численного моделирования детонационных волн.
Моделируемая система представляет собой замкнутую идеально теплоизолированную емкость постоянного объема, в начальном состоянии заполненную смесью газообразных веществ, способных вступать в систему суммарно-экзотермических реакций. Состав смеси, ее начальная температура и давление задаются. В такой системе при любых начальных условиях выделяющееся в ходе реакций тепло приводит к прогрессирующему разогреву смеси, закапчивающемуся тепловым взрывом через Длительность периода индукции зависит как от начальных условий и теплофизических свойств среды, так и от параметров химической реакции (константы скорости, энергии активации, теплового эффекта). В ходе процесса меняется химический состав реагирующей смеси и ее температура.
Задача об адиабатическом тепловом взрыве впервые была рептетта О.М. Тодесом в 1933 г. [1]. В пределе простой одностадийной сильно экзотермической реакции выражение для оценки выглядит так [2]
= £ М. *. ,,,
Здесь с, р - теплоемкость и плотность реагирующей смеси; К- универсальная газовая постоянная; То - начальная температура реагирую щей смеси; Q,E — тепловой эффект и энергия активации реакции; го - т.н. предэкспонент, имеющий размерность обратную времени и зависящий от начальной концентрации реагентов и порядка реакции [2, 3].
Применение этой простой формулы весьма ограничено, поскольку большинство процессов горения и взрыва (в частности - в газовой фазе) имеют сложные кинетические механизмы разветвленных цепных реакций. Они включают в себя комплексы из десятков, сотен и даже тысяч химических реакций (пг) с участием десятков-сотен реагентов (пс), демонстрируют кинетическое поведение, принципиально отличающееся от поведения простых реакций [4] и, чаще всего, не могут быть описаны единой простой брутто-реакцией, пригодной для моделирования процесса горения в широком диапазоне изменения концентраций и температур. Для математического моделирования цепных реакций приходится численно решать
пс
дающие высокой степенью математической жесткости [о, 6]. Аналитические решения здесь невозможны, и численное моделирование является единственно возможным инструментом для исследования подобных задач.
1. Постановка задачи
Объект моделирования - изохорттый адиабатический реактор представляет собой замкнутую, идеально теплоизолированную емкость постоянного объема V, в которой находится химически реагирующая смесь переменного состава из пс газообразных компонентов
— исходных веществ, промежуточных и конечных продуктов реакций, а так же инертных веществ в реакциях не участвующих. В таком реакторе могут одновременно протекать пг элементарных химических реакций, механизм и кинетические константы для которых різ-
В0СТНЫ.
В начальный момент времени і = 0 задаются: начальная температура То, давление Ро: состав смеси задается набором начальных мольных долей компонентов:
С
N
і = 1
■ Пс
(2)
х
где Ог, N - мольная концентрация ^Ог = и количество молей г-го компонента сме-
си. Обозначим как Оз = ^^ Ог - суммарную концентрацию и N3 = N - суммар-
гг
ттое количество молей в смеси. Соответственно начальные значения будем обозначать как Х0г, Оог, N0г, Ооз, N03
Важным допущением является отсутствие начальных градиентов концентраций, температуры и давления в объеме реактора - т.е. в начальный момент во всех точках реактора они одинаковы. Поскольку система идеально изолирована, па стенках реактора градиенты концентраций, температуры и давления равны нулю. Следствием этого является отсутствие градиентов на протяжении всего процесса - концентрации, температура и давление будут меняться синхронно во всех точках реагирующей смеси. Математическая модель реактора может быть описана системой обыкновенных дифференциальных уравнений с начальными условиями, в которой единственной независимой переменной будет время t.
2. Вывод основных соотношений для идеального газа
Используем общепринятые принципы построения моделей макрокиттетики и химической газодинамики [2-5].
а) Уравнение баланса массы. В замкнутой системе постоянного объема масса т и
плотность р вещества не меняются
т =^ т =^ тог = то; р = ^ Р“ =^ Рог = Ро; Р = у; 1 = 1 •••Пс• (3)
гг гг
При этом, массы отдельных компонентов тг и их плотности рг могут меняться во времени из-за протекания химических реакций, по соотношения (3) должны выполняться. Добавим еще несколько соотношений, связывающих массы и плотности с концентрациями
Рг = ЦгОг; р = ^ Цг Ог; т = V ^ ЦгОг, (4)
гг
где Цг ~ молярная масса г-го компонента смеси.
Начальные масса и плотность смеси определяются по начальным мольным долям, температуре и давлению с помощью уравнений состояния. Подробнее этот момент будет обсуждаться ниже.
б) Уравнение баланса массы г-го компонента смеси. Как уже отмечалось выше, массы отдельных компонентов тг могут меняться во времени из-за протекания химических превращений. В замкнутом безградиептпом реакторе это единственный фактор, влияющий па баланс массы, поэтому скорость изменения тг может быть получена из соотношения
7 ПГ
-тг .
"-Г = £ ^, (5)
3=1
где Ггз - скорость производства или расходования массы г-го компонента в реакции. Запишем ^-го реакцию в виде обобщенного стехиометрического уравнения
у13 А1 + у2 3А2 + ••• + уПс,3АПс ^ у1,3А1 + у2,3А2 + ••• + уПс,3АПс (6)
Здесь Аг - химический символ г-го реагента (например Н2, СО2, СН4, СН2О и т.д.); у,3, у,,3 -стехиометрические коэффициенты г-го вещества в реакции, которые показывают сколько молей (или молекул) данного реагента вступают в реакцию как исходные вещества (,г,3) или получаются как продукты реакции ^3Стехиометрические коэффициенты веществ, не участвутопщх в реакции равны 0.
Скорость газофазной реакции определяется законом действующих масс (зависимость от концентраций) и законом Аррениуса (зависимость от температуры) [2, 3] и может быть выражена как
К = г,Т*' ехр(-3 П . (7)
' ' г=1
Здесь Х3, Ь3, Е3 - кинетические константы для расчета скорости реакции.
Скорость производства и расходования массы г-го компонента в j-vt реакции выражается
гг,3 = К3 (уг,3 — уг3 Цгу, (8)
а скорость изменения массы г-го вещества во всем комплексе реакций получается суммиро-
ванием вкладов каждой реакции в производство или потребление данного вещества
7 Пг
= „,гУ^Ш3 (,/,3 - Уг3 • О)
3=1
С учетом (4) данное выражение можно записать в более удобном для практического применения виде
dC Пг
-; {=1 ■■■Пс (10)
j=i
в) Уравнение баланса энергии. Поскольку энергетический обмен реактора с окружающей средой запрещен, суммарная энергия рассматриваемой системы не меняется, что может быть сформулировано в виде постоянства суммарной энтальпии. Но обмен энергией между компонентами системы происходит непрерывно, в частности в экзотермических химических реакциях химическая энергия переходит в тепловую, в эндотермических - наоборот.
Уравнение баланса энергии запишем как
Пс
^Cihi = const, (11)
i=1
здесь hi - удельная мольная энтальпия г-го вещества). Уравнение сохранения энергии удобнее записать в дифференциальной форме
d ± Cihi = t (Ci+ hif ) =0. (12)
i=l i=l 4 7
В этом выражении dCi/dt определено выше зависимостью (10), а производную энтальпии преобразуем как
dhi dhi dT dT
dt dT dt Cpi dt ’
здесь Cpi - изобарная молярная теплоемкость г-го вещества; (12) принимает вид
dT dC dt + i dt
E(CiCpi '§ + h =0. (13)
i=1
В случае идеального газа энтальпия и теплоемкость зависят только от температуры и из (13) получаем уравнение для контроля текущей температуры газовой смеси.
Пс Г АС-
Т - £ Ь Т ^
* = ------------------1' (14)
(т)]
г=1
При высоких давлениях необходимо применять уравнение состояния реального газа, а энтальпия и теплоемкость будут зависеть и от давления. Особенности построения энергии для этого случая будут обсуждаться ниже.
Зависимости Нг (Т) и сРг (Т)для всех веществ задаются в виде термодинамических полиномов. Зависимости для Нг (Т) получены аппроксимацией табличных зависимостей энтальпии из [7] полиномами 3-го порядка, работающими в диапазоне температур от 300 до 4000 К. Зависимости для сРг (Т) получены дифференцированием по температуре полиномов для энтальпии.
г) Задание начальных концентраций. Для решения поставленной задачи - определения периода индукции данная система должна быть дополнена начальными условиями, задающими начальные концентрации и температуру
Ь = 0, Т = То, Сг = С0г. (15)
В предыдущем разделе было отмечено, что начальное состояние реактора определяется иначе - задается начальная температура, давление и мольные доли реагентов (2). Переход от этой системы данных к (15) возможен с использованием уравнения сохранения массы (4) и применяемого уравнения состояния газовой смеси.
В случае идеального газа різ уравнения состояния
Ро = ^ ЯТа = Сов ЯТа (16)
выражаем начальную суммарную концентрацию
Сов = ЯТ-. (17)
и начальные концентрации всех компонентов
Соі = ХіСов; І = 1 ...Пс. (18)
Таким образом, мы можем решить сформулированную выше задачу численно. Математическая модель адиабатического реактора описывается уравнением баланса энергии (14) для определения температуры реагирующей смеси Т и пс уравнениями (10) для мольных концентраций каждого різ реагентов. На каждом птаге рептетшя мы получаем вектор теку-Сі Т
можтго рассчрітать по формуле (16)
P = RT^TCi (19)
i=1
tad
тый в теоррш теплового взрыва кррітеррій - момент достріжєпрія максрімума скорострі тепло/ cT
выделения I max ——
\ dt
3. Выбор уравнения состояния реального газа
Для корректного опрісатірія процессов в газовой фазе пррі вьісокріх давлетшях тіеобхо-дрімо прріметіять уравпетше состоятірія, учрітьівающее срільі взаимодействрія молекул пррі ріх сблріжетірірі, в частпострі срільі отталкріватгая. Уравпетше должно отвечать следутоїтцім тре-боватіРіям:
а) Уравпепріе должно быть работоспособным пррі очень вьісокріх давлетшях - как мрітірі-мум до 2000 МПа.
б) Структура этого уравпетгая должна позволять ему гтібко прріспосаблріваться к сложной смєсрі газов переменного состава, что важно пррі опрісатгарі хрімрічєскрі реагрірутощей среды.
в) Это же уравпепріе должно быть пррігодпо pi для отдельного опрісатірія любого рітідрі-врідуальпого компонента смєсрі.
Исходя різ этртх требоватіРій, різ мтіогообразрія уравпетіРій состоятірія, прріметіяемьіх для опрісатірія реального газа [8, 9], для дальнейшего пррімепетшя было выбрано уравпепріе Беккера - Крістяковсого - Врільсотіа (общепррітіятое сокращетше BKW), [10J
PVm xK вхк
m = 1 +-------x-----е (T+°)а (20)
rt + (T + в)а к }
Здесь Ут - молярный объем газовой смеси (величина обратная суммарной концентрации С^ коэффициент К для рассмотренной выше смеси реагентов определяется выражением
Пс
К = ^ Сгкг, (21)
г=1
где кг - т.н. коволюмы индивидуальных компонентов смеси; Сг - их концентрации.
Как показано в [8, 10] существует несколько методик определения значений коволтомов (геометрический, через потенциал Леппарда - Джонса, экспериментальный). В [10] приведена таблица значений коволтомов, определенных разными методами - разброс величин для одного вещества достигает десятков процентов. В дальнейших расчетах использовались значения, приведенные в табл. 1. Большая часть значений взята из [10], предпочтение отдавалось экспериментально определенным значениям, при отсутствии таковых брались наименее противоречивые значения, определенные по потенциалу Леппарда - Джонса или геометрические коволюмы. Для ряда веществ данные отсутствуют, и значения коволтомов были рассчитаны самостоятельно по геометрической методике, описанной в [10].
Таблица 1
Коволюмы индивидуальных веществ кг, см3/моль
В-во кі В-во кі В-во кі В-во кі
Н2 120 Н С 400 Н2 О2 590 С2Н3 1150
02 300 С2Н6 1100 СН 472 С2Н4 1100
N2 380 Н 60 СН2 520 С2Н5 1100
Н2О 250 О 120 СНз 520 СНО 860
СО 380 ОН 410 С2Н 950 СН2О 815
2 о С 600 2 О Н 590 С2 Н2 1013 СН3О 860
Эмпирические параметры а, в, X, 0 подбираются по экспериментальным данным, в частности по адиабатам Гтоготтио для продуктов детонации. В [10] приведена таблица рекомендуемых параметров уравнения ВК^¥ для продуктов детонации различных взрывчатых веществ. Все эти наборы параметров, а также некоторые другие, встречающиеся в литературе по детонации, были протестированы тта задаче моделирования зависимости молярных объемов отдельных газов от давления и температуры (см. ниже). Сравнение полученных зависимостей с экспериментальными данными из книги Циклиса [11] показало, что наиболее удовлетворительный результат дал единственный набор параметров а = 0, 25, в = 0, 3. X = 1, 0 = 0, рекомендованный Хэллфордом [8].
Уравнение состояния в этом случае принимает более простой вид
РУт К 03 К
К— = 1 + ^ е ^ • (22)
Это уравнение, называемое так же уравнением Хэлфорда - Кистяковсого - Вильсотта, рекомендовано в [8] для решения задач детонации при давлениях до 20 000 МПа.
Данные соображения и привели к однозначному выбору уравнения (22) в качестве объекта моделирования.
Оценка работоспособности уравнения проводилась по следующей методике. Если записать уравнение (22) тте для смеси газов, а для объема заполненного чистым веществом г-го сорта, мы получим уравнение состояния индивидуального газа
РУті Сікі 0,3сікі кі °’3к*95 , ,
—— = 1+-------— е Т0,25 = 1+-------і--е ушіТ°’26 ( 23)
КТ Т° > 25 е УтіТ° > 25 е • { ]
Его можно преобразовать к виду стандартного алгебраического уравнения, трансцендентного относительно мольного объема Vmi.
PV ■ к °’3к^
PVmi - 1-------k^-5еVmiT02 = 0. (24)
RT VmT02 К J
Численно рептая это уравнение при различных сочетаниях Р и Т можно получить таблицу функции Vmi = f (Р, T, ki) и сравнить полученные результаты с экспериментальными данными, приводимыми в литературе. Так в [11] приведен набор экспериментальных зависимостей мольных объемов для ряда газов. Для решаемой ниже задачи моделирования горения метан - воздушной смеси есть данные только для пяти компонентов - N2, H2, O2, CO2, CH4. Решение уравнения (24) производилось с помощью встроенной в пакет Mathcad функции root и с помощью программного модуля, реализующего быстро сходящийся метод секущих-хорд, который далее был использован в основной Фортран - программе.
Для водорода рачетпые данные хорошо согласуются с литературными данными; в остальных случаях уравнение (24) дает несколько завышенные значения мольных объемов, причем с повышением давления расхождение данных уменьшается.
4. Особенности модели адиабатического реактора в случае реального газа
4.1. Расчет начального состава смеси в реакторе с помощью уравнения состояния реального газа
Вернемся к задаче расчета начальных концентраций при заданном давлении Ро, температуре То и составе смеси в виде мольных долей реагентов xoi- Порядок расчета таков.
а) Рассчитываем мольные объемы для всех компонентов смеси, для которых xoi = 0. решая трансцендентное уравнение (24) при Ро, То.
б) Далее руководствуемся следующими соображениями. Если общее (заранее неизвестное) количество молей Nos поместить в реактор объемом V, то количество молей г-го
вещества будет равно Noi = XoiNos. С другой стороны, это количество молей будет
при заданных условиях занимать объем Vi = NoiVmi, а сумма этих объемов должна быть равна объему реактора
Пс Пс Пс
V 'У ' NoiVmi ^ ^ XoiNoS Vmi = NoS^~''y XoiVmi.
i=l i=l i=l
Пс
Поделив правую и левую часть этого соотношения на V, имеем Cos xoi Vmi — 1;
i=1
что при уже известных мольных объемах дает выражение для расчета суммарной концентрации
Cos = -n—1---------------------------------------. (25)
і=1
в) Определяем начальные концентрации всех компонентов Соі = ХіЄов; і = 1 •• -пс-
Очевидно, что из-за худшей сжимаемости реального газа, при росте давления эти концентрации будут меньше, чем при том же давлении и температуре в идеальном газе. Это
должно привести, в соответствии с формулой (7), к замедлению реакций и увеличению периода индукции, причем этот эффект будет усиливаться с ростом давления.
Если по ходу решения необходимо рассчитывать текущее давление для вычисления поправок к энтальпии и теплоемкости, то для этого используется уравнение (22)
4.2. Влияние давления на энтальпию и теплоемкость реального газа
Как показано в [12], уравнение, выражающее зависимость энтальпии от давления при постоянной температуре выглядит следующим образом
Разделяя переменные и интегрируя (27) от Ро до Р получаем выражение для изменения энтальпии за счет роста давления (в нашем случае Ро = 0,1 МПа)
Здесь, как и раньше, Утг - мольный объем г-го вещества. Поскольку Утг определяется решением трансцендентного уравнения (24), аналитическое вычисление интеграла и производной в правой части (29) невозможно. Их можно вычислить с помощью численных методов. При отработке методики с помощью математического пакета МаШсаА, были использованы встроенные функции численного дифференцирования и интегрирования.
Полученные по (29) значения оказались выше, чем литературные оценки [11], полученные с помощью аналогичного уравнения, по с использованием уравнения состояния Тэйта и методов графического дифференцирования. При этом общий характер зависимостей практически одинаков.
Поправку к теплоемкости г-го вещества за счет высокого давления можно получить путем численного дифференцирования выражения (29) по температуре
Существенным недостатком численного расчета поправок (29), (30) является необходимость многократного расчета мольных объемов Уті, которые в свою очередь, получаются итерационным решением трансцендентного уравнения (24). Поэтому применение высокоточных встроенных методов численного дифференцирования и интегрирования приводит к неоправданно большому времени расчета, которое ощутимо даже па высокопроизводительных компьютерах. С целыо снижения вычислительных затрат было произведено исследование возможности применения простых вычислительных схем мепыттей точности, по требующих минимума вычислений. После серии численных экспериментов выбор был оста-
НОВЛ6Н Неї СЛвДуЮЩИХ МбТОДеІХ.
і=1
(26)
(27)
н(Р,Т) - н(Ро,т) = р у - т(д-^ ар.
(28)
Переписав это соотношение для парциальных мольных энтальпий, получаем
АЫ(Р,Т) = Ні(Р, Т) - Ні(Ро,Т)= Р Уті - т^^УтА ар; і = 1
■;р0 І V /р.
і = 1 •••Пс• (29)
І = 1 • ••Пс^
(30)
а) Для вычисления производных в формулах (29), (30) с сохранением минимум 6 верных значащих цифр хоропто подходит центрально-разностная формула 2-го порядка точности при вариации температуры па 0,1К.
б) Для численного вычисления интеграла в (29) оптимальным оказалось применение метода паивысптей алгебраической точности Гаусса с 8 узлами, обеспечивающего при минимуме вычислений 4-5 точных значащих цифр, что для данной задачи вполне достаточно.
Именно эти методы и были применены в основной Фортран-программе.
4.3. Модификация уравнения баланса энергии
Для случая реального газа уравнения (10) остаются без изменения, а вывод уравнения баланса энергии должен быть пересмотрен. В формуле (12) появляются дополнительные слагаемые, учитывающие влияние давления. Производную энтальпии преобразуем так
dhi
dt
dhi dT dhi dP dT dt + dP dt
dT . dhi (dP dT dP dCi dP \dT dt + dCi dt
Уравнение энергии принимает вид
Пс
Е
i=1
C , c + dhi dP\ dT + (h + C dhi dP\ dCi i 1 Cpi + dP dT dt M i + i dP dC\ dt
(31)
(32)
В случае реального газа уравнение для контроля текущей температуры газовой смеси будет выглядеть тэж
dhi dP \ dCi ' dP dCi/ dt
dT
dt
Е
i=1
hi (T) +Ahi(P,T) + Ci
Ci
i=1
* (T) + AcPi(P,T) + § §
(33)
Апалрітріческріе выражетшя для расчета производных давления по температуре pi концентрации получаются дифференцированием (26)
dP K
dT = Cs R 11 + T°>25
1 -0,“U+Ti)
К = J2 Ciki.
(34)
i=1
. + C ki ( 0 3K
1 + CsК I To25
І = 1 ...Ur.
(35)
Значение для производной энтальпии по давлению рассчитываем численно в соответствии с приведенным выше выражением (27)
= Vmi - T dP ) T mi
dVm,
~dT
(36)
P
5. Реализация основного вычислительного алгоритма
Для решения всех задач, описанных ниже, использована система программирования Compaq Digital Visual FORTRAN v. 6.6. При программировании модулей для расчета правых частей уравнений химической кинетики pi матррщы Якобрі Ріспользовалась тєхпологрія прямых безрщдекспых ВЬІЧРІСЛеіІРІЙ, ЧТО существенно повышает скорость ВЬІЧРІСЛеіІРІЙ. Исходные тексты соответствутопщх модулей па ФОРТРАНЕ геперрірутотся автоматріческрі с помощью спецріальпо разработанной программы.
0.3K
Решение жесткой системы дифференциальных уравнений выполнялось с использованием солвера ЭТИТ, основанного па многошаговом методе Гира переменного порядка, специально адаптированного для решения задач химической кинетики [6]. Отт имеет весьма эффективный блок адаптивного выбора птага. Позволяет использовать для вычисления матрицы Якоби внешнюю пользовательскую процедуру, хотя имеет и встроенные средства для ее численного расчета. Длительный опыт применения данной программы показал высокую эффективность и устойчивость при решении широкого круга задач теории горения.
Для расчетов использовалась широко известная и многократно апробированная кинетическая схема окисления легких углеводородов, предложенная в [13]. Использован несколько укороченный вариант (первые 78 обратимых реакций). Исключены ненужные для решаемой задачи реакции окисления метанола. Поскольку обратимую реакцию удобно рассматривать как две независимые односторонние реакции, принятая схема включает 156 реакций, 23 активных реагента (Н2, 02, ЩО, СО СО2, СН4, С2Н3, И, О, ОН И02, Н202, СИ, СН2, СН3. С2Н, С2Н2, С2Н3, С2Н4 С2Н5, СНО, СН2О, СН3О) и один инерт N2.
Таким образом, решается система різ 25 дріфферетщріальттьіх уравпетгай (24 концентрации + температура), пс = 24; пг = 156.
Для реїттепрія задачрі расчета адріабатріческого перргода ріпдукцрірі прріметіялрісь 3 модрі-фрікацрірі программы:
Вариант 1 реалрізует модель процессов в идеальном газе па основе уравпетшй (10), (14);
Вариант 2 реалрізует модель процессов в реальном газе па основе уравпетшй (10), (33);
Вариант 3 реалрізует комбинированную модель процессов в газе тіа основе уравпетшй рідє-альттого газа (10), (14), по с учетом начального сжатрія в соотвєтстврірі с уравпетшем (25) для реального газа; влріяттріє давлетгая тіа эптальпрпо рі теплоемкость тіе учріттлва-ется.
В связрі с тем, что что варріатітьі 1 рі 3 программы факттіческрі отлрічатотся только способом определеттрія начального состава, бьістродействріе этртх программ почтрі тіе отлрічается. Варріатіт 2 программы, ттапротрів, пррі каждом обращеттрш к модулю расчета правых частей дріффереттцріальттьіх уравпетгай, выполняет массу дополтгательттых оператцій, связанных с расчетом текущего давлетгая; поправок тіа давлетгае к эптальпршм 24 веществ; чріслєтітіьім рітітегрріроватіріем рі дріфферетщріроватіріем; расчетом поправок в уравпепрш эттергрга. Пррі этом прріходрітся многократно перевьічріслять мольный объем каждого компонента смєсрі рітерацріотітіьім рептетіРіем трансцендентного уравпетгая (24). В результате, этот варріатіт программы решает задачу тіа порядок медленнее остальных.
6. Результаты расчетов
Рассчрітьіватотся адріабатріческріе перргоды ріпдукцрірі для разных составов метаті-кріслородтіьіх рі метаті-воздупттіьіх смесей пррі заданном начальном давлепрш. Расчеты про-водрілрісь для следутощріх начальных давлетгай: 0,1, 1, 10, 100, 1000 МПа. Пррі заданном давлетіРіРі, расчеты проводятся одтгам потоком последовательно для серрга начальных температур от 800 до 2500 К через 100 К.
а) Сравнение периодов индукции в идеальном и реальном газе для стехиометрической метан-кислородной смеси (2О2 + СН4)
Использовалрісь варріатітьі 1 рі 2 программ. Результаты расчетов прріведетіьі тіа рріс. 1 рі, выборочно, в табл. 2. Хорошо врідтто, что пррі тірізкріх давлетгаях перргоды ріпдукцрірі в рідє-альттом рі реальном газах отлрічатотся тіезпачрітельтіо. Одтіако, пачрітіая с давлетгая 10 МПа,
разница становится достаточно ощутимой. Начиная с 100 МПа, она видна даже на графиках в логарифмических координатах. Для давлений порядка 1000 МПа период ріпдукцрірі в реальном газе в несколько раз болытте, чем в Рідеальтіом. Пррі росте начальной температуры отпосрітельпая разтшца в перргодах ріпдукцрірі Рідеальтіого рі реального газов уменьшается.
3 -Ю 4 4-Ю 4 5 -Ю 4 6 10 4 7 -Ю 4 8 -Ю 4 9 -Ю 4 0.001 0.0011 0.0012
1/Т, к1
Рис. 1. Стехиометрическая метан-кислородная смесь. Сравнение результатов расчетов для идеального (сплошная линия) и реального газов (пунктир)
б) Сравнение периодов индукции в идеальном и реальном газе для стехиометрической метан - воздушной смеси (2О2 + СН4 + 7, 52^)
Расчеты проводрілрісь для тех же начальных давлетгай тємрі же варріантамрі программы. Результаты расчетов прріведетіьі па рріс. 2 рі, выборочно, в табл. 3. Здесь тепдетітція та же самая - пррі тірізкріх давлетшях перргоды ріпдукцрірі отлрічатотся тіезпачрітельпо, а тіа-чрітіая с с давлетгая 10 МПа, разтшца в перргодах ріпдукцрірі в Рідеальтіом рі реальном газах статіоврітся достаточно ощутрімой. Пррі этом отпосрітельпая разтшца пррі одріттаковой температуре несколько болытте, чем в случае метап-кріслородпой смєсрі. Пррі росте начальной температуры отпосрітельпая разтшца в перргодах ріпдукцрірі Рідеальтіого рі реального газов также уменьшается.
в) Сравнение периодов индукции в двух моделях реального газа для стехиометрической метан - кислородной смеси (2О2 + СН4) и стехиометрической метан
- воздушной смеси (2О2 + СН4 + 7, 52^)
С целью выделетгая домрпшрутопщх факторов, влріятоіцріх тіа перргод ріпдукцрірі в реальном газе, была выполнена серрія атталогрічттьіх расчетов варріатітом 3 программы, которая
Таблица 2
Результаты расчета периода индукции, с для смеси 202 + СН4
Р, МПа Модель газа Температура, К
800 900 1000 1200 1500 2000 2500
0,1 Идеальный Реальный Комбин. 202,0 202,2 202,2 7,945 7.951 7.952 0,577 0,577 0,577 0,013 0,013 0,013 4,00е-4 3,87е-4 4,00е-4 1,42е-5 1,42е-5 1,42е-5 2,39е-6 2,39е-6 2,39е-6
1 Идеальный Реальный Комбин. 22,17 22,38 22,37 0,819 0,830 0,826 0,056 0,057 0,057 1,20е-3 1,21е-3 1,21е-3 4,28е-5 4,ЗЗе-5 4,30е-5 1,66е-6 1,67е-6 1,67е-6 2,47е-7 2,48е-7 2,48е-7
10 Идеальный Реальный Комбин. 2,50 2,72 2,71 0,089 0,097 0,095 5,86е-3 6,28е-3 6,25е-3 1,08е-4 1,09е-4 1,14е-4 3,60е-6 3,73е-6 3,77е-6 1,85е-7 1,91е-7 1,91е-7 2,74е-8 2,80е-8 2,80е-8
100 Идеальный Реальный Комбин. 0,2659 0,4478 0,4413 9,254е-3 0,0151 0,0147 6,03е-4 9,47е-4 9,26е-4 1,08е-5 1,59е-5 1,56е-5 3,23е-7 4,45е-7 4,40е-7 1,77е-8 2,28е-8 2,25е-8 3,12е-9 3,78е-9 3,72е-9
1000 Идеальный Реальный Комбин. 0,0271 0,1312 0,1264 9,43е-4 4,255е-3 4,050е-3 6,11е-5 2,60е-4 2,46е-4 1,08е-6 4,09е-6 3,87е-6 3,20е-8 1,04е-7 1,00е-7 1,75е-9 4,79е-9 4,61е-9 3,27е-10 7,63е-10 7,ЗЗе-10
Таблица 3
Результаты расчета периода индукции, с для смеси 202 + СН4 + 7, 52К2
Р, МПа Модель газа Температура, К
800 900 1000 1200 1500 2000 2500
0,1 Идеальный Реальный Комбин. 831,1 832.0 832.0 34,80 34.84 34.84 2,722 2.724 2.724 0,068 0,068 0,068 1,78е-3 1,79е-3 1,79е-3 5,30е-5 5,30е-5 5,30е-5 8,54е-6 8,54е-6 8,55е-6
1 Идеальный Реальный Комбин. 88,18 89,07 89,06 3,482 3,511 3,514 0,258 0,260 0,260 6,61е-3 6,65е-3 6,65е-3 2,35е-4 2,36е-4 2,36е-4 6,65е-6 6,69е-6 6,67е-6 9ДЗе-7 9,22е-7 9Д5е-7
10 Идеальный Реальный Комбин. 10,07 11,03 11,01 0,376 0,410 0,407 0,026 0,028 0,028 5,43е-4 5,52е-4 5,77е-4 2,12е-5 2,24е-5 2,23е-5 9,11е-7 9,40е-7 9,36е-7 1Д1е-7 1Д4е-7 1ДЗе-7
100 Идеальный Реальный Комбин. 1,099 1,891 1,872 0,040 0,066 0,065 2,69е-3 4,31е-3 4,24е-3 5,12е-5 7,80е-5 7,60е-5 1,64е-6 2,38е-6 2,38е-6 1,05е-7 1,39е-7 1,ЗЗе-7 1,58е-8 1,92е-8 1,87е-8
1000 Идеальный Реальный Комбин. 0,113 0,571 0,554 4,11е-3 1,95е-2 1,85е-2 2,73е-4 1,23е-3 1,16е-3 5,12е-6 2,10е-5 1,93е-5 1,56е-7 5,73е-7 5Д6е-7 8,31е-9 2,52е-8 2,83е-8 1,72е-9 4,64е-9 4,37е-9
Рис. 2. Стехиометрическая метан воздушная смесь. Сравнение результатов расчетов для идеального (сплошная линия) и реального газов (пунктир)
фактически является моделью идеального газа, за исключением того, что начальные концентрации ораделятотся на основе уравнения реального газа. В таблицах 2 и 3 показаны результаты сравнения этой серии расчетов (в графе «комбинированная модель:») с зависимостями, полученными в пунктах а) и б). Результат оказался довольно неожиданным - из таблиц ясно, что максимальные расхождения результатов, полученных вариантами 2 и 3 программы, не превышают нескольких процентов, да и то только при высоких давлениях.
7. Анализ полученных результатов
а) Очень слабое влияние довольно больших по величине поправок к энтальпии и теплоемкости, показанное в предыдущем разделе, па первый взгляд выглядит весьма неожиданным. Именно энтальпии определяют тепловые эффекты реакций, и их изменение, казалось бы, должно приводить к изменению температуры в соответствии с уравнениями (14), (33). Температура же очень сильно влияет на скорость реакции в соответствии с законом Аррениуса (7), что может вызвать изменение периода индукции. Объяснить наблюдаемое отсутствие такого сильного эффекта можно тем, что тепловой эффект реакции определяется, в соответствии с законом Гесса [3], разностью энтальпий продуктов реакции и ее исходных веществ. В случае, если при изменении давления энтальпии веществ меняются синхронно, разность энтальпий будет сохраняться.
Одттако вопрос о синхронности изменения энтальпий остается открытым. В соответствии с формулой (29) поправка па давление рассчитывается через мольные объемы компонентов, а они все же меняются по-разному, т.к. у реагентов коволтомы отличаются по величине. Детальное исследование промежуточных результатов вычислений показало, что дополнительные слагаемые в числителе и знаменателе правой части уравнения (33) имеют тот же порядок, что и основные слагаемые в уравнении (14). Одттако, определяющие скорость изменения температуры соотношения числителя и знаменателя, вычисляемые по формулам (14) и (33) сохраняют близкие по величине значения в ходе всего периода индукции. Поправки взаимно компенсируют друг друга.
Этот эффект дает весьма полезный результат. Если пас устраивает определение периода индукции с ошибкой в несколько процентов, да и то только при очень высоких давлениях, можно в расчетах не вычислять очень трудоемкие поправки на давление для термодинамических параметров.
б) В свете сказанного выше, легко выделить доминирующий фактор, определяющий влияние давления тта период индукции. Это начальное сжатие реагирующее смеси. Запишем уравнение состояния в обобщенном виде
РУ = г(Ро, Тс),
МЯТ
где г - коэффициент сжимаемости. Для идеального газа г = 1, а для реального газа при высоком давлении г > 1.
Выразим начальную суммарную концентрацию при заданных Ро, То
°03 = ^У~ = г (Ро, То)яТо ■ (37)
Отсюда становится поттяттто, что начальная концентрация в идеальном газе (при прочих равных условиях) всегда будет больше, чем в реальном газе при высоком давлении, причем с повышением давления эта разница будет возрастать. Увеличение температуры несколько уменьшает г (см., например, уравнение (22)) и повышает начальную концентрацию.
Уменьшение начальной концентрации в реальном газе, рассчитанной по формуле (25), эквивалентной (37), приведет к пропорциональному уменьшению концентраций промежуточных продуктов и, в соответствии с законом действующих масс, к замедлению скорости всех реакций. В следствие этого период индукции увеличивается. Иметттто это и определяет все особенности поведения периодов индукции в реальном и идеальном газе, отмеченные в предыдущем разделе, в том числе и уменьшение относительного отклонения периодов индукции с ростом температуры.
в) Не следует забывать, что рассматриваемая модель адиабатического реактора является весьма идеализированной. При низких начальных температурах и давлениях периоды индукции составляют соттти и тысячи секунд, и обеспечить тта практике действительную адиабатичттость реактора можтто только с некоторой степенью приближенности.
Напротив, при высоких температурах и давлениях периоды индукции, становятся чрезвычайно малыми. На практике такие начальные условия создаются ударным сжатием газа, которое выполняется тте мгттоветттто. Если период индукции становится соизмерим с временем сжатая, или даже становится меттытте его, данная модель уже тте будет работать, т.к. тта ход реакций будет влиять скорость ттарастания температуры и давления.
Работа проводилась при финансовой поддержке проекта 2012 066 - Г327 ЮУрГУ.
Литература
1. Тодес, О.М. «Адиабатический» тепловой взрыв / О.М. Тодес // ЖФХ. - 1933. - Т.4. № 1. - С. 71-78. (Сб. Теория горения и взрыва. - М.: Наука, 1981. С. 195-202).
2. Франк-Камепецкий, Д.А. Диффузия и теплопередача в химической кинетике / Д.А. Франк-Камепецкий. - М.: Наука, 1987. - 502 с.
3. Математическая теория горения и взрыва / Я.Б. Зельдович, Г.И. Барепблатт, Б.Б. Либ-рович, Г.М. Махвиладзе. - М.: Наука, 1980. - 478 с.
4. Варпатц, К). Горение. Физические и химические аспекты, моделирование, эксперименты, образование загрязняющих веществ / К). Варпатц, У. Маас, Р. Диббл. - М.: ФИЗ-МАТЛИТ, 2003. - 352 с.
5. Оратт, Э. Численное моделирование реагирующих потоков / Э. Оратт, Дж. Борис. - М.: Мир, 1990. - 660 с.
6. Полак, Л.С. Вычислительные методы в химической кинетике / Л.С. Полак, М.Я. Голь-деттберг, А.А. Левицкий. - М.: Наука, 1984. - 280 с.
7. Термодинамические свойства индивидуальных веществ. В 4 т. /Л.В. Гурвич и др. - М.: Наука, 1978-1982. - 4 т.
8. Гирптфельдер, Дж. Молекулярная теория газов и жидкостей / Дж. Гирптфельдер,
Ч. Кертис, Р. Берд. - М.: Изд-во иностр. лит., 1961. - 930 с.
9. Рид, Р. Свойства газов и жидкостей / Р. Рид, Дж. Праусттиц, Т. Шервуд. - Л.: Химия, 1982. - 592 с.
10. Мейдер, Ч. Численное моделирование детонации / Ч. Мейдер. - М.: Мир, 1985. - 384 с.
11. Циклис, Д.С. Плотные газы / Д.С. Циклис. - М.: Химия, 1977. - 168 с.
12. Калашников, Я.А. Физическая химия веществ при высоких давлениях / Я.А. Калашников. - М.: Выспт. птк., 1987. - 241 с.
13. Вестбрук, Ч. Применение химической кинетики для определения критических параметров газовой детонации / Ч. Вестбрук, П. Уртьев /'/ ФГВ. - 1983. - Т. 19, №6. - С. 65-76.
Валерий Константинович Рябинин, кандидат физико-математических паук, доцент, кафедра «Вычислительная механика сплошных сред», Южно-Уральский государственный университет (г. Челябинск, Российская Федерация), rowan@math.susu.ac.ru.
Юрий Михайлович Ковалев, доктор физико-математических паук, профессор, кафедра «» ситет (г. Челябинск, Российская Федерация), yuin_kov@inail.ru.
Bulletin of the South Ural State University. Series «Mathematical Modelling, Programming & Computer Software:»,
2013, vol. 6, no. 1, pp. 56-71.
MSC 80A25, 80A30, 97M10
Mathematical Modelling of Adiabatic Induction Period for Methane-Oxygen Mixtures in a Wide Range of Initial Pressure and Temperature
V.K. Ryabinin, South Ural State University, Chelyabinsk, Russian Federation,
rowan@math.susu.ac.ru,
Yu.M. Kovalev, South Ural State University, Chelyabinsk, Russian Federation,
yum _ ко v@mail. ru
The subject of this work is the features of thermal explosion modeling in an isochoric adiabatic reactor for wide range of initial pressures and temperatures, with reference to the maximal kinetic mechanism of a branched - chain reaction for methane - oxygen and methane - air mixtures. Two mathematical models are constructed on the basis of gase law equations for ideal and actual gases. The amendments to an enthalpy and heat capacity for components of a gas mixture due to high-pressure in real gas are taken into account. Results of calculations for these two classes of models shown, that at high pressure the induction period in real gas is longer, than in ideal at the same starting conditions. Set, that the dominant factor, defining this effect, is difference in density of real and ideal gas after initial compression.
Keywords: mathematical modelling, combustion and explosion, chain reactions, high pressure.
References
1. Todes O.M. «Adiabaticheskiy> teplovoy vzryv [«Adiabatio Heat Explosion]. ZhFKh [J. of Physical Chemistry], 1933, vol. 4, no. 1, pp. 71-78.
2. Frank-Kamenetskii D.A. Diffusion and Heat Transfer in Chemical Kinetics. Plenum Press, 1969.
3. Zel’dovich Ya.B., Barenblatt G.I., Librovich B.B., Makhviladze G.M. Matematicheskaya teoriya goreniya i vzryva [The Mathematical Theory of Combustion and Explosion]. Moscow, Nauka, 1980.
4. Warnatz J., Maas U., Dibble R.W. Combustion. Physical and Chemical Fundamentals, Modeling and Simulations, Experiments, Pollutant Formation. Springer, 2001.
5. Oran E.S., Boris ,J.P. Numerical Simulation of Reactive Flow. Elservier Science Publishing Co., Inc., New York, 1987.
6. Polak L.S., GoPdenberg M.Ya., Levit.zkiy A.A. Vychislitelnye metody v himicheskoy kinetike [Computational Methods in Chemical Kinetics]. Moscow, Nauka, 1984.
7. Gurvich L.V. Termodinamicheskie svoistva individual 'nykh veshchestv. In 4 vo^s [Thermodynamic Properties of Individual Substances]. Moscow, Nauka, 1978-1982.
8. Hirschfelder J.O., Curtiss Ch.F., Bird R.B. Molecular Theory of Gases and Liquids. John Wiley and Sons, Inc., Ney York, Chopinan and Hall, Lirri., London, 1954.
9. Reid R.C., Prausnitz J.M., Sherwood Т.К. The Properties of Gases and Liquids. McGraw-Hill, Inc., 1977.
10. Mader Ch.L. Numerical Modeling of Detonations. University of California Press, 1979.
11. Tsiklis D.C. Plotnye gasy [Dense Gases]. Moscow, Chimiya, 1977.
12. Kalashnikov Ya.A. Fizicheskaya khimiya veshchestv pri vysokikh davleniyah [Physical Chemistry of Substances Under High Pressure]. Moscow, Vysshaya shcola, 1987.
13. Westbrook Ch.K., Urtiew P.A. Use of Chemical Kinetics to Predict Critical Parameters of Gaseous Detonations. Combustion, Explosion and Shock Waves, 1983, vol. 19, no. 6, pp. 753-766.
Поступила в редакцию 7 ноября 2012 г.
