Химическая физика, 2020, T. 39, № 12, стр. 29-36

ВВЕДЕНИЕ

Изучение воспламенения, горения и детонации в газовых смесях обусловлено широким распространением этих процессов в различных технологических, промышленных, энергетических установках, а также в двигательных системах. Проведение масштабных экспериментов сопряжено с огромными материальными затратами, поэтому на многих этапах исследований активно применяется численное моделирование [1]. Данная расчетная процедура позволяет выявить определяющие особенности протекания различных физико-химических процессов и удачно прогнозировать экспериментальные исследования. В ряде случаев с помощью численного моделирования удается прогнозировать оптимальные условия проведения экспериментов [2]. Особенно это важно для изучения сложных пограничных процессов, когда наблюдаются переходы очагового воспламенения в распространение пламени, а затем уже – горения в детонацию [3‒6].

Благодаря быстрому развитию возможностей вычислительной техники и совершенствованию математического аппарата для численного моделирования с высокой точностью динамики газа появляется возможность более глубокого анализа и понимания процессов воспламенения, горения и детонации. Это приводит к значительному увеличению размеров детального кинетического механизма (ДКМ) по количеству компонент и реакций. В работах [7, 8] представлена динамика увеличения ДКМ за последние 30 лет. На основании ее анализа делается вывод о том, что наблюдается резкое увеличение не только числа реакций и компонент, составляющих различные ДКМ, но и количество самих ДКМ, разрабатываемых различными научными группами. В итоге ДКМ содержат сотни компонент и несколько тысяч элементарных реакций. Причин такого структурного увеличения наполняемости ДКМ и роста количества самих ДКМ – три. Первая связана с тем, что появляются новые данные об элементарных химических реакциях. Большинство этих количественных данных в настоящее время получаются с помощью квантовомеханических расчетов [9, 10]. Вторая причина обусловлена методикой конструирования ДКМ, когда блоки новых реакций добавляются к уже существующим ДКМ [11, 12], тем самым расширяя последние. Третья относится к вопросам мониторинга процессов горения и пиролиза, когда развитие методик регистрации позволяет получать более точную кинетическую информацию о концентрациях все большего числа атомов, радикалов и молекул, включая электронно-возбужденные состояния некоторых из них [13‒15]. Мониторинг в совокупности с численным моделированием горения дает возможность значительно улучшить энергоэффективность и экологическую безопасность сжигания различных углеводородных топлив.

Значительное расширение ДКМ по числу компонент и реакций приводит к проблеме их прямого использования в газодинамических расчетах с учетом химической кинетики, обременяя расчеты увеличением времени счета и накоплением вычислительных ошибок [16‒18]. В этой связи приходится проводить оптимизацию ДКМ с целью уменьшения количества компонент и реакций, получая в итоге редуцированный кинетический механизм (РКМ) [7, 8, 18, 19] и/или глобальный кинетический механизм [20, 21]. Методики редуцирования ДКМ в основном сводятся к уменьшению числа реакций [16‒18]. В то же время ранг матрицы Якоби определяется количеством компонент в ДКМ [1]. Если расчет скоростей реакций сводится к арифметическим действиям с матрицами, то интегрирование проводится по количеству компонент. В этой связи разумно понижать количество компонент, а вместе с ними и реакций, в которых они участвуют. В настоящее время математический аппарат редуцирования по числу реакций достаточно хорошо представлен различными методами, описанными в работах [16‒18]. Однако методик редуцирования числа компонент в литературе представлено значительно меньше. В основном они сводятся к трем методам.

1. Первым из них является метод квазистационарных состояний [17]. Этим методом в работе [22] было проведено редуцирование ДКМ, состоящего из 28 компонент и 114 реакций, до РКМ, содержащего 25 реакций и 15 компонент для стехиометрической смеси С₂Н₂/О₂ в диапазонах температуры 900 < T < 2000 K и давления 0.01 < T < 5 МПа. Данный метод требует от исследователя хорошего знания химической кинетики конкретного процесса горения или пиролиза. Он удачно применяется в тех случаях, когда число компонент в ДКМ незначительное.

2. Для ДКМ, в которых число компонент может измеряться сотнями, применяется метод проб и ошибок (trial-and-error method) по анализу отклика к заданному физическому параметру [1, 16, 18].

3. В последнее время наибольшее распространение получили различные варианты метода анализа графа прямых связей (DRG), основанного на определении влияния отдельных компонент на группу компонент, определяющих поведение параметров химических превращений [18, 19, 23, 24].

Во всех трех методах важным является выбор физикохимических параметров, по которым проводится редуцирование, и тех интервалов времени, в которых проводится анализ системы обыкновенных дифференциальных уравнений (ОДУ), определяющих химическую кинетику. Как следствие, редуцирование по методу DRG становится отдельной научно-исследовательской задачей, что делает этот процесс трудноосуществимым в автоматическом режиме [18, 24].

В настоящей работе апробирована методика автоматического редуцирования числа компонент и реакций ДКМ по задержке воспламенения τ для стехиометрических смесей ацетилена и метана с кислородом. Выбор объектов исследования обусловлен необходимостью продемонстрировать универсальность предлагаемого подхода. Ацетилен обладает самыми высокими энергетическим потенциалом и детонационной способностью среди других углеводородов, в то время как метан является его антиподом [25‒28].

МЕТОДИКА РЕДУЦИРОВАНИЯ ДКМ

Сложности расчетов с использованием ДКМ напрямую связаны с решением прямой кинетической задачи, описываемой системой ОДУ, которая в общем случае является нелинейной [18, 29, 30]. Данная система ОДУ определена законом действующих масс, согласно которому изменение во времени концентрации i-го компонента X_i равно

где $R_{j}^{ + }$ и $R_{j}^{ - }$ – скорости образования и расходования компонента X_i в j-той реакции; $N_{R}^{i}$ – количество реакций, в которых участвует X_i; $R_{j}^{ + }$ и $R_{j}^{ - }$ определяются как здесь k_j есть константа скорости j-той реакции; ν – порядок реакции; l – индекс, соответствующий реагенту в левой части j-той реакции; ξ – число расходующихся в j-той реакции компонентов X_i. Обычно порядок реакции ν ≤ 3, и, соответственно, ξ ≤ 2. В случае гомогенного горения в замкнутом пространстве при адиабатических условиях или в системах с фиксированной скоростью тепловых потерь система ОДУ (1) дополняется дифференциальным уравнением для скорости изменения температуры: где C_p – теплоемкость системы при постоянном давлении, N_R – общее число реакций, R_j – скорость j-той реакции, χ – коэффициент теплопередачи из системы во внешнюю среду, S и V – площадь поверхности и объем системы, T – текущая температура системы, T₀ – температура окружающей среды, ${{\Delta }_{r}}H_{j}^{0}$ – стандартная молярная энтальпия j‑той реакции. Таким образом, в правой части уравнения (1) имеет место функция Ψ, в общем случае зависящая во времени от Т, Р, [X_i], k_j, $~{{\Delta }_{r}}H_{j}^{0}$, C_p и ряда других параметров.

Система ОДУ (1) в совокупности с уравнением (4) принадлежит к классу так называемых “жестких” систем, которые характеризуются большим разбросом собственных значений матрицы Якоби [29], определяемой функцией Ψ. Последнее приводит к измельчению шага интегрирования и в итоге – к накоплению вычислительной ошибки и увеличению времени счета [1, 18].

Для уменьшения числа компонент в данной работе использовался метод проб и ошибок, а для уменьшения числа реакций – методика оценки вклада каждой реакции в изменение для каждого X_i правой части системы (1). Такой подход является вариантом анализа скоростей образования (rate-of-production analysis) [18]. На каждом шаге интегрирования, согласно (2) и (3), определялись вклады $R_{j}^{ + }$ и $R_{j}^{ - }$ в скорость образования, $\sum\nolimits_j^{N_{R}^{i}} {R_{j}^{ + }} ,$ и расходования, $\sum\nolimits_j^{N_{R}^{i}} {R_{j}^{ - }} ,$ i‑той компоненты по всем реакциям, участвующим в этих процессах. Далее по критериям [31]

определялись те реакции, которые предстояло исключать из ДКМ. Здесь ε – параметр малости реакции, который принимался всегда равным 0.001. В случае медленных процессов, когда наблюдается квазистационарность по компонентам, критерии (5) объединяются с критерием где ${{R}_{j}}$ = $R_{j}^{ + }$ ‒ $R_{j}^{ - }$.

Процедура редуцирования состояла в следующем. Строилась сетка значений Т_l и Р_l с фиксированным количеством точек l, покрывающих весь диапазон значений температуры и давления. По значениям этой сетки проводился расчет τ(T_l, P_l) для ДКМ. Затем путем последовательного исключения каждого i-го компонента определялось значение τ_‒i(T_l, P_l), которое сравнивалось с τ(T_l, P_l) по критерию

где δ – допустимое относительное отклонение. Если значение τ_‒i(T_l, P_l) удовлетворяло критерию (7) для всех l, то i-тый компонент удалялся из ДКМ вместе со всеми реакциями, в которых он участвовал. После прохождения всех точек сетки Т_l и Р_l проводилось редуцирование оставшихся реакций согласно критериям (5) и (6). Из рассмотрения по критерию (7) следует сразу исключить исходные реагенты и те активные атомы, радикалы и молекулы, взаимодействие которых определяет цепной разветвленный процесс самовоспламенения углеводородов. К ним относятся активные атомы О и Н, радикалы ОН, НО₂, НСО, молекулы Н₂, СО, СО₂, СН₂О, Н₂О₂, Н₂О [1, 17, 25].

В соответствии с критериями (5)‒(7) была написана программа на языке FORTRAN 77, соединенная с программным модулем CHEMKIN III [32]. Все вычисления в данной работе проводились в предположении V = const. Процедура редуцирования осуществлялась в диапазонах температуры 1000 < T < 2100 К для смеси 21.43% C₂H₂ + 53.575% O₂ в аргоне (φ = 1) и 1300 < T < 2000 К для смеси 10% CH₄ + 20% O₂ в аргоне (φ = 1). Для смеси C₂H₂/O₂ расчеты проводились для трех давлений: 0.02, 0.1 и 1.1 МПа. Для смеси CH₄/O₂ значение давления составило 0.65 МПа.

АНАЛИЗ РАСЧЕТОВ ПО РАЗЛИЧНЫМ ДКМ

Прежде чем приступить к редуцированию, мы провели обзор ДКМ, представленных в литературе. Их сравнение осуществлялось на основе описания температурной зависимости задержки воспламенения τ. Выбор τ в качестве параметра обусловлен его высокими обобщенными характеристиками для анализа процесса самовоспламенения углеводородов [25]. Поскольку программный модуль CHEMKIN III требует использования термодинамических параметров для всех содержащихся в ДКМ компонент в формате NASA [33], то численное моделирование проводилось только с представленными в литературе ДКМ [11, 12, 34‒39], содержащими соответствующие термодинамические файлы данных. Под задержкой самовоспламенения горючих смесей в расчетах принималось значение τ, соответствующее максимальному росту давления [25].

На рис. 1 представлены значения τ, полученные с использованием ДКМ из работ [11, 35‒39]. Из этого рисунка видно, что при Р = 0.02 и 0.1 МПа ДКМ из [11, 36, 38] предсказывают близкие значения τ на всем диапазоне температур, и только при Р = 1.1 МПа ДКМ из работы [36] предсказывает большее значение τ при T < 1200 K. В то же время значения τ, рассчитанные по ДКМ из [35], дают самые короткие, а по ДКМ из [37] – наболее длинные задержки воспламенения при всех значениях Р и Т. Значение τ из [39] меньше величин τ из работ [11, 36, 38], но демонстрирует аналогичную с ними температурную зависимость при всех трех значениях Р. Как показывает рис. 1, для Р = 0.02 и 0.1 МПа зависимость τ(Т) практически не изменяется для ДКМ из [11, 36, 38, 39], в то время как для ДКМ из [37] наблюдается меньшая зависимость τ(Т) при температуре выше 1400 К. При Р = 1.1 МПа уже для ДКМ из [11, 36, 38, 39] проявляется такое же уменьшение τ(Т), что и для ДКМ из [37]. В то же время значения τ, рассчитанные по ДКМ из [35], при всех трех значениях Р имеют одинаковую зависимость от температуры, слабо уменьшающуюся с ростом Т. Такое поведение τ(Т) приводит к тому, что если при T = = 2100 К значения τ уменьшаются с ростом Р, то при T = 1000 К для всех упомянутых ДКМ, кроме ДКМ из [35], наблюдается либо слабая, либо обратная зависимость τ(Р). Подобные изменения зависимостей τ(Р) от температуры, рассчитанные по ДКМ из работ [11, 36‒39], были получены в работе [22] для смесей ацетилена с кислородом в том же диапазоне температуры при давлениях от 0.01 до 5 МПа. Следует отметить, что с понижением температуры и ростом давления растет и расхождение значений τ, рассчитанных по разным ДКМ (рис. 1).

В целом, рис. 1 показывает, что расчеты в рамках различных ДКМ дают различающиеся результаты: как по характеру, так и в абсолютных значениях температурной зависимости задержки воспламенения τ для смеси ацетилена с кислородом. При этом различия в абсолютных значениях τ могут превышать один порядок величины для разных ДКМ. Выбор ДКМ, предназначенного для редуцирования, был основан на следующих соображениях. Согласно рис. 1, ДКМ из [11, 36, 38, 39] показывают аналогичные зависимости τ(Т) и τ(Р) в промежутке между крайними значениями τ, рассчитанными по ДКМ из [35, 37]. Расчеты значений τ по ДКМ из [11, 38] предсказывают близкие результаты для всех трех давлений во всем температурном диапазоне. Из этих двух механизмов для редуцирования был выбран ДКМ из [38], содержащий 57 компонент и 405 реакций.

РЕДУЦИРОВАНИЕ ДКМ И ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Редуцирование ДКМ из [38] для смеси C₂H₂−O₂−Ar с использованием сетки из девяти точек при давлениях 0.02, 0.1 и 1.1 МПа и температурах 1000, 1500 и 2100 К осуществлялось при значении δ = 20% и свелось к РКМ(25;124), состоящему из 25 компонент и 124 реакций. Сопоставление результатов расчетов задержек воспламенения ДКМ из [38] и РКМ(25;124) представлено в табл. 1, где значения τ демонстрируют хорошее совпадение расчетов по ДКМ [38] и РКМ(25;124) на всей сетке значений Т и Р.

Расчеты τ, полученные с помощью РКМ(25;124), сопоставлялись с экспериментальными данными, представленными в работах [40‒43] (рис. 2 и 3). Рисунок 2 демонстрирует хорошее соответствие описания температурной зависимости задержки воспламенения τ, рассчитанной по РКМ(25;124), экспериментальным данным из работ [40, 41]. Здесь же представлены и результаты расчетов τ, проведенные с помощью ДКМ из работ [11, 12, 34, 39]. Хотя на рис. 2 представлены задержки воспламенения для сильно разбавленной аргоном смеси C₂H₂/O₂, тем не менее видно, что расчет по ДКМ из [11] близок к расчетам по РКМ(25;124), как и в случае условий для рис. 1. Близки к ним и расчеты по ДКМ из работ [12, 39]. Расчет по ДКМ из [34] показал большее значение τ во всем температурном диапазоне.

На рис. 3 представлено сравнение рассчитанных по РКМ(25;124) и ДКМ из работ [11, 12, 39] и экспериментально наблюдаемых температурных зависимостей задержки воспламенения стехиометрических смесей C₂H₂/O₂/Ar по данным из работ [42, 43]. Поскольку при близких составах смесей и величинах давления значения τ, полученные в работах [42] и [43], хорошо согласуются между собой, расчеты проводились для условий, приведенных в [42]. Как и в случае рис. 2, видно, что РКМ(25;124) достаточно хорошо описывает экспериментальные результаты во всем температурном диапазоне. Детальные КМ из работ [11, 39] также показали хорошее описание температурной зависимости τ при T > 970 K. В то же время с понижением температуры эти ДКМ дали несколько большее значение τ, чем в [42]. Детальный КМ из [12] дает хорошее описание этой зависимости только при T > 1100 К. С понижением температуры значение τ по ДКМ из [12] увеличивается и становится в три раза больше, чем в [42].

Апробация методики редуцирования ДКМ из работы [38] применительно к смеси 10% CH₄ + + 20% O₂, разбавленной аргоном, была осуществлена с δ = 30% при Р = 0.65 МПа в температурном диапазоне 1300 < T < 2000 K. Использовались три значения температуры: 1300, 1600 и 2000 К. В результате был получен РКМ(16; 60), состоящий из 16 компонент и 60 реакций. Таблица 2 демонстрирует хорошее соответствие расчетов по РКМ(16; 60) и ДКМ из [38]. Видно, что для смеси CH₄/O₂/Ar наибольшее расхождение расчетов по РКМ и исходному ДКМ из [38] наблюдается при Т = 2000 К. Сравнение расчетов по РКМ(16;60) и ДКМ из работ [11, 12, 34, 39] представлено на рис. 4 при описании экспериментально наблюдаемых значений τ из [31]. Расчеты проводились для φ = 1, поскольку в [31] было показано, что изменение φ от 0.33 до 3 в указаном диапазоне температуры слабо влияет на задержку воспламенения, изменяя ее приблизительно в 2 раза (рис. 4). Такая слабая зависимость τ(Т) от φ позволяет применительно к углеводородовоздушным смесям проводить редуцирование ДКМ при φ = 1. Из рис. 4 видно, что расчеты задержек воспламенения τ по РКМ(16; 60) и ДКМ из работ [11, 12, 34, 39] достаточно хорошо согласуются с данными из [31].

Варьируя величину δ, можно изменять степень редуцирования. Так, при изменении значения δ с 20% на 10% для смеси С₂Н₂/О₂ ДКМ из работы [38] свелся к РКМ, состоящему из 37 компонент и 235 реакций. Прямой зависимости между количеством компонент и величиной δ не существует, поскольку каждый компонент вносит с блоком своих реакций вклад в формирование значения τ(Р, Т) не только по абсолютной величине Δτ, но и по знаку +/‒. Аналогичным влиянием обладают и группы взаимосвязанных компонентов, что демонстрируют различные методы DRG [18]. Однако общая тенденция уменьшения числа компонент с увеличением δ наблюдается.

В представленной работе редуцирование ДКМ осуществлялось по параметру τ, что было обусловлено моделированием процесса самовоспламенения ацетилена и метана в смесях с кислородом. В других задачах в качестве параметров редуцирования могут выступать иные характеристики процессов физико-химических превращений. Такими параметрами (η) могут служить нормальная скорость распространения пламени, профили концентраций исходных реагентов, продуктов горения и пиролиза, суммарное тепловыделение [25]. В случае моделирования излучения и ионизации в реагирующих смесях целесообразно принимать в качестве η как сами электронно-возбужденные состояния частиц и зондовые токи, так и концентрации атомов, радикалов и молекул, определяющие эти процессы [13‒15, 40, 44]. Предложенный метод редуцирования ДКМ можно применять и в условиях, когда нарушается термодинамическое равновесие, например для сверхзвуковых течений, где требуется учет влияния неравновесной колебательной и химической кинетики и необходимо оценивать ее влияние на параметры потока [45]. Для этого нужно дополнить систему ОДУ (1) соответствующими уравнениями поуровневой кинетики [45‒47]. Аналогичным образом можно расширить область применения предложенного метода редуцирования ДКМ для гетерогенных физико-химических процессов с учетом каталитических реакций [19, 30, 48].

ЗАКЛЮЧЕНИЕ

Проведенное численное моделирование стехиометрической смеси С₂Н₂/О₂ в широких диапазонах температуры и давления позволило выбрать из представленных в литературе оптимальный ДКМ для дальнейшего его редуцирования. Определив τ в качестве основного характеристического параметра процесса самовоспламенения углеводородных топлив, с помощью метода проб и ошибок удалось эффективно средуцировать исходный ДКМ применительно к смесям С₂Н₂/О₂ и СН₄/О. Разработанный программный модуль позволил в автоматическом режиме по параметру δ(τ) получить РКМ для каждой из исследуемых смесей углеводородов. Апробация полученных РКМ по представленным в литературе экспериментальным данным показала их хорошую валидацию по описанию задержек воспламенения. Развитие методики редуцирования числа компонент ДКМ требует разработки сложного математического аппарата для поиска минимумов функционала по уравнениям (1) и (4) в автоматическом режиме. Успех решения проблемы редуцирования ДКМ во многом будет определяться выбором параметра η и балансом точности описания процесса с определением значения δ. Представленная работа показала, что и без исследования функции Ψ можно эффективно проводить редуцирование ДКМ.

Работа выполнена в рамках Программы фундаментальных исследований Российской академии наук на 2013–2020 гг. по теме ФИЦ ХФ РАН № 49.23 (номер государственной регистрации ЦИТИС: АААА-А18-118031590088-8).