Химическая физика, 2019, T. 38, № 8, стр. 69-79

ВВЕДЕНИЕ

Прямое численное моделирование турбулентного горения предполагает неэмпирический подход на основе применения детального кинетического механизма (ДКМ) и учета поля пульсационных скоростей. Первые работы по математическому моделированию турбулентного горения появились достаточно давно [1, 2], однако общепринятых методов математического описания физико-химических процессов в турбулентном пламени до сих пор в литературе нет (см. [3–9]). Существующие подходы, как правило, включают всякого рода “замыкающие” гипотезы, основанные на экспериментальных наблюдениях. Поскольку область применимости той или иной гипотезы ограничена, такие полуэмпирические подходы неуниверсальны.

Сегодня наиболее перспективным универсальным подходом к неэмпирическому теоретическому описанию турбулентного горения считают прямое численное моделирование (ПЧМ, в англоязычной литературе – Direct Numerical Simulation (DNS)), при котором в рассмотрение включены все основные особенности трехмерного турбулентного реагирующего течения с полным спектром турбулентных пульсаций скорости, полным набором исходных, промежуточных и конечных химических компонентов с их индивидуальными термохимическими свойствами и свойствами молекулярного переноса, а также с адекватными граничными условиями [7–9]. Численное интегрирование определяющих уравнений течения проводят на расчетных сетках, обеспечивающих пространственное разрешение турбулентных вихрей колмогоровского масштаба, используя схемы высокого порядка аппроксимации.

Несмотря на значительные успехи в развитии ПЧМ турбулентного горения, на этом пути есть еще много проблем. Так, при сравнении самых “свежих” приведенных в литературе решений для скорости распространения турбулентного пламени, полученных с помощью ПЧМ (см., например, работу [10]), с известными экспериментальными данными выясняется, что их точность все еще заметно уступает точности, достигаемой при решении задач ламинарного горения (авторы работы [10] получили величины скорости распространения турбулентного пламени в метановоздушной смеси, в 3–5 раз большие экспериментальных значений, описанных в литературе). В теоретических работах в качестве причин расхождения результатов расчетов и экспериментов часто называют неучтенные эффекты, вызванные неопределенными граничными условиями в реальных экспериментах. Кроме того, применение ПЧМ к реальным турбулентным пламенам требует очень больших вычислительных ресурсов, в связи с чем в расчетах все еще приходится использовать различные упрощения, влияющие на точность решения.

В данной работе используется альтернативный подход к ПЧМ турбулентного горения гомогенной газовой смеси, который ранее был предложен в работе [11] в двумерном приближении. Вместо численного решения всех уравнений, определяющих распространение турбулентного пламени в реагирующем газе, в [11] предложено решать только уравнения переноса скалярных величин – концентраций реагентов и энергии – в искусственном (“синтетическом”) поле турбулентности, характеризуемом заданной (постоянной) среднеквадратичной интенсивностью пульсаций скорости с заданными (постоянными) интегральными пространственными и временными масштабами. При этом предполагалось, что химический процесс описывается одной глобальной реакцией, а само распространение пламени не влияет на характеристики синтетического поля турбулентности в предпламенной зоне.

В отличие от [11], в работе [12] химический процесс описывался уже не одной глобальной реакцией, а ДКМ с участием активных центров. В [12] результаты по моделированию турбулентного горения в двумерном приближении были получены для одного конкретного состава водородно-воздушной смеси в простейших условиях синтетического поля стационарной, однородной и изотропной турбулентности. В работе [13] методика из [12] была усовершенствована и применена для расчета распространения турбулентного пламени в водородно-воздушных гомогенных смесях разного состава и при разных давлениях.

Цель данной работы – переход от моделирования турбулентного горения в двумерном приближении к физически адекватной трехмерной задаче.

МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ

Постановка задачи в общем виде сформулирована в работе [13], но здесь целесообразно ее повторить. Система уравнений, определяющих распространение турбулентного пламени, основана на уравнениях Навье–Стокса в трехмерной постановке, уравнении сохранения энергии и уравнениях неразрывности для всех химических компонентов в смеси идеальных газов [14]:

где $t$ – время; $\rho $ – плотность; ${\mathbf{v}}$ – вектор скорости; $p$ – статическое давление; $e$ – внутренняя энергия; ${\mathbf{q}}$ – вектор молекулярного потока тепла; ${\mathbf{P}}$ – тензор сил давления, ${{{\mathbf{f}}}_{i}}$ – вектор силы тяжести, действующей на единицу массы i-го вещества; ${{Y}_{i}},$ ${{{\mathbf{V}}}_{i}},$ ${{h}_{i}}$ и ${{\omega }_{i}}$ – массовая доля, вектор скорости диффузии, удельная энтальпия и скорость химического превращения i-го вещества соответственно; $h_{i}^{0}$ – стандартная энтальпия образования i-го вещества; $N$ – число веществ в реагирующем газе; ${{R}^{0}}$ – универсальная газовая постоянная, ${{\mu }_{i}}$ – молекулярный вес i-го вещества; $T$ – температура; ${{c}_{{p,i}}}$ – теплоемкость i-го вещества при постоянном давлении; $\nabla $ – дифференциальный оператор; начальные условия в свежей смеси отмечены индексом “0”.

Предполагается, что система (1) дополнена ДКМ окисления горючего, термохимическими данными каждого вещества ($h_{i}^{0},$ ${{c}_{{p,i}}}$) и соответствующими соотношениями для ${{{\mathbf{f}}}_{i}},$ ${\mathbf{q}},$ ${\mathbf{P}},$ ${{{\mathbf{V}}}_{i}}$ и ${{\omega }_{i}},$ а также начальными и граничными условиями. В результате решения задачи должны быть получены скорость стационарного распространения ${{u}_{t}}$ (скорость турбулентного горения) и структура турбулентного пламени в многокомпонентном реагирующем газе. Чтобы упростить решение задачи, введем следующие основные допущения:

1. Область течения имеет простейшую геометрию; турбулентность стационарная, однородная и изотропная.

2. Давление постоянно ($p = {{p}_{0}}$), а влиянием силы тяжести можно пренебречь. Благодаря этим допущениям задача существенно упрощается, так как отпадает необходимость решать уравнение сохранения количества движения.

3. Тепловой поток ${\mathbf{q}}$ определяется только молекулярной теплопроводностью (лучистым теплообменом можно пренебречь).

4. Эффекты термодиффузии малы.

5. Реагирующая смесь сильно разбавлена инертным газом (азотом), так что диффузионные потоки всех химических компонентов определяются законом Фика с коэффициентом бинарной диффузии.

С учетом допущений 2–5 система дифференциальных уравнений (1) сводится к виду

где $\lambda $ – коэффициент молекулярной теплопроводности газа, а ${{D}_{i}}$ – коэффициент бинарной диффузии i-го вещества в азоте. Отметим, что благодаря второму предположению плотность $\rho $ в уравнении (2) – функция только температуры и состава смеси.

Для дальнейшего упрощения задачи предположим, что решение системы (2) стремится к некоторому стационарному решению, для которого

Поскольку нас интересует именно такое стационарное решение, уравнение (3) можно использовать вместо уравнений неразрывности смеси, т.е. вместо первого уравнения в системе (2). В этом случае постановка задачи допускает дальнейшее упрощение. Вектор мгновенной местной скорости ${\mathbf{v}} = (u,v,w)$ в (3), где $u,\,\,{v},\,\,w$ – компоненты вектора, можно представить в виде суммы вектора средней скорости ${\mathbf{V}} = (U,V,W)$ (здесь прописные буквы означают среднее значение) и вектора пульсаций скорости ${\mathbf{v}}{\text{'}} = (u{\text{'}},{v}{\text{'}},w{\text{'}})$ (штрихи относятся к пульсациям скорости): ${\mathbf{v}} = {\mathbf{V}} + {\mathbf{v}}{\text{'}}.$

Тогда (3) можно преобразовать к виду

Основываясь на первом допущении (об однородной и изотропной турбулентности), предположим далее, что

и, следовательно,

В области течения с простейшей геометрией (см. допущение 1) всегда существует усредненное направление распространения турбулентного пламени, и можно устремить одну ось пространственных координат в этом направлении. Тогда выражение (6) можно проинтегрировать и получить

где ${\mathbf{B}}$ – постоянный поток вещества. С учетом (7) систему (2) можно привести к виду

Последние члены в обоих уравнениях системы (8) определяют перенос вещества и энергии с турбулентными пульсациями скорости. Если вектор ${\mathbf{B}}$ известен, то производные по времени в этих уравнениях можно положить равными нулю, чтобы получить структуру стационарного турбулентного пламени. Однако, как правило, значение вектора ${\mathbf{B}}$ заранее неизвестно. Поэтому, чтобы получить его значение, необходимо решить нестационарную систему (8), описывающую движение волны горения до достижения некоторой постоянной скорости. В этом случае конвективные члены в системе (8) можно опустить, а система уравнений примет вид

Наконец, можно упростить и последние члены в двух уравнениях (9), используя (5) и предполагая, что

где ${{\nabla }_{x}} = {\partial \mathord{\left/ {\vphantom {\partial {\partial x}}} \right. \kern-0em} {\partial x}},$ ${{\nabla }_{y}} = {\partial \mathord{\left/ {\vphantom {\partial {\partial y}}} \right. \kern-0em} {\partial y}}$ и ${{\nabla }_{z}} = {\partial \mathord{\left/ {\vphantom {\partial {\partial z}}} \right. \kern-0em} {\partial z}}.$ Интегрирование уравнений (10) дает

После подстановки соотношений (11) в систему (9) окончательно получим

Таким образом, как и в работах [11–13], вместо численного решения системы (1) предлагается решать систему (12), состоящую только из уравнений переноса скалярных величин – концентраций всех N реагентов и энергии – на заданном синтетическом поле турбулентности, характеризуемом среднеквадратичной интенсивностью пульсаций скорости v' с заданными интегральными пространственными и временны́ми масштабами. При трехмерном моделировании в системе координат (x, y, z) уравнения (12) в развернутом виде и с учетом калорического уравнения состояния будут выглядеть следующим образом:

Система (13) состоит из N + 2 уравнений для N + 2 переменных ($N$ веществ с массовой долей ${{Y}_{i}},$ температурой $T$ и плотностью $\rho $). Для замыкания системы используем следующие соотношения для ${{c}_{p}},$ $\lambda $ и ${{D}_{i}}{\text{:}}$

где a₁–a₅ – коэффициенты полиномов, а индекс “in” относится к инертному разбавителю (азот).

Для скорости химических реакций используется соотношение [14]

где $\nu _{{i,k}}^{ + }$ и $\nu _{{i,k}}^{ - }$ – стехиометрические коэффициенты для i-тых веществ, являющихся реагентами и продуктами k-й реакции, соответственно; ${{A}_{k}}$ – предэкспоненциальный множитель k-й реакции; ${{\alpha }_{k}}$ – показатель степени, определяющий температурную зависимость предэкспоненциального множителя k-й реакции; ${{E}_{k}}$ – энергия активации k-й реакции; $M$ – общее количество химических реакций и ${{X}_{j}}$ – мольная доля j-го вещества.

Уравнения (13) содержат параметры стационарной, однородной и изотропной турбулентности реагирующего газа (три компоненты вектора пульсационной скорости: $u_{0}^{'},\,\,{v}_{0}^{'},\,\,w_{0}^{'}$). Как и в работах [11–13], компоненты вектора пульсационной скорости получаем путем статистического разыгрывания методом Монте-Карло, предположив, что мгновенные местные компоненты вектора пульсационной скорости удовлетворяют нормальному распределению Гаусса φ, а вихревая структура турбулентности описывается экспоненциально затухающими пространственной (R') и временнóй (R) корреляционными функциями:

Соотношение (15) для краткости записано только для x-составляющей мгновенной пульсационной скорости u' со среднеквадратичным значением u, причем индекс “0” в обозначениях опущен. Такое же соотношение используется для разыгрывания y- и z-составляющих мгновенной пульсационной скорости ${v}'$ и w' с соответствующими среднеквадратичными значениями ${v}$ и w. В (15) и (16) $\sigma $ – среднеквадратичное отклонение пульсаций скорости по направлению x, y или z; r_k – расстояние в физическом пространстве; L' и τ – пространственный и временнóй масштабы турбулентности. В соответствии с первым допущением (см. выше) о стационарности и однородности поля турбулентности, масштабы τ и L' приняты постоянными.

Среднеквадратичная длина пространственного вектора пульсационной скорости $\bar {u}$ определяется по среднеквадратичным значениям его компонент u, ${v}$ и w соответственно по осям x, y и z с помощью соотношения

Подчеркнем, что задание вектора пульсационной скорости $\bar {u}$ и пространственного масштаба турбулентности L' равносильно заданию диссипации турбулентной энергии ε ~ ${{\bar {u}}^{3}}$/L' и колмогоровского масштаба турбулентности l₀ ~ a^3/4ε^–1/4, где a = λ/ρc_p – коэффициент температуропроводности, который для газов близок к коэффициенту кинематической молекулярной вязкости.

Пример статистического разыгрывания поля пульсационных скоростей u' по одной из координатных осей (по оси х) методом Монте-Карло представлен на рис. 1. Получающиеся при этом зависимости среднеквадратичных компонент скорости u, ${v}$ и w от времени приведены на рис. 2. Как видно, компоненты скорости u и ${v}$ строго соответствуют задаваемым исходным значениям, а компонента w, за счет которой выдерживается закон неразрывности, в самом начале процесса оказывается больше задаваемого значения, однако затем приближается к u и ${v}.$ Поэтому для выполнения условия изотропности турбулентности самое начало процесса (соответствующее области до вертикальной штриховой линии на рис. 1) в проведенных расчетах исключалось из рассмотрения. Граничные и начальные условия для системы (13) рассмотрены ниже для конкретной расчетной области.

МЕТОДИКА РАСЧЕТОВ

Расчетная область представляет собой параллелепипед, верхняя и нижняя грани которого – плоские квадраты. В начальный момент времени t = t₀= 0 объем содержит два слоя, разделенных плоской поверхностью, параллельной основанию и отстоящей от него на некоторое расстояние. Нижний слой содержит продукты горения, верхний – свежую смесь. Нижняя и верхняя грани удалены от разделительной поверхности на достаточно большое расстояние, чтобы обеспечить постоянство давления в системе за все время расчета. Боковые грани – непроницаемые стенки со скольжением потока. Пламя распространяется снизу вверх. На рис. 3 для примера показана трехмерная поверхность распространяющегося пламени, искривленная турбулентными пульсациями скорости.

На рис. 4 представлена двумерная схема расчетной области в выбранной системе координат (x, y) при некотором значении z. Начальное положение поверхности пламени в момент времени t = t₀= 0 показано горизонтальной штриховой линией. При t > 0 положение пламени условно представляется плоскостью, параллельной основанию параллелепипеда, выше которой находится только свежая смесь, а ниже – только продукты реакции. Подчеркнем, что это – условный фронт пламени, поскольку сама зона реакции в пламени имеет конечную толщину. Таким образом, для определения скорости распространения пламени условно принимается, что свежая смесь занимает пространство с температурой Т < T_m, а продукты реакции – пространство с температурой T > T_m. Значение T_m равно среднему арифметическому между начальной температурой свежей смеси, T₀, и термодинамически равновесной температурой продуктов сгорания, T₁:

На рис. 4 мгновенное положение фронта пламени в момент времени t = t_i > 0 – след истинной поверхности пламени на плоскости ху – показан тонкой сплошной линией, а положение условного фронта пламени – горизонтальной штриховой линией.

Задачу о распространении турбулентного пламени решали с помощью численного интегрирования системы (13) с дополнительными соотношениями методом переменных направлений [15, 16]. Для решения уравнений по каждому из направлений (x, y, z) применялась неявная разностная схема с равномерным шагом по пространству и линеаризацией нелинейных “источниковых” слагаемых на верхнем слое. Эта схема имеет первый порядок точности по времени и пространству. Расчетная область пространства (x, y, z) имела размер 1 × 0.5 × 1 см³. Все расчетные ячейки имели форму куба размером 0.01 × 0.01 × 0.01 см³ или 0.005 × 0.005 × 0.005 см³. Полное количество расчетных ячеек составило 500 000 и 4 000 000 соответственно. Шаг интегрирования по времени изменялся в зависимости от числа итераций, но не превышал 10^–6 с.

Важная операция – моделирование синтетического поля турбулентности $v_{0}^{'} = (u_{0}^{'},{v}_{0}^{'},w_{0}^{'})$ в уравнениях (13). Моделирование полей пульсационных скоростей $u{\text{'}}$(x, y, z, t), $v{\text{'}}$(x, y, z, t) и w'(x, y, z, t) проводилось путем разыгрывания их возможных значений при заданных параметрах турбулентности u, ${v},$ w, L' и τ. Процедура разыгрывания предусматривала выполнение закона неразрывности. Соответствующая методика описана в Приложении.

Ниже представлены используемые здесь начальные и граничные условия. Начальные условия –

Граничные условия –

РЕЗУЛЬТАТЫ РАСЧЕТОВ

Как и в [13], в данной работе рассматривается турбулентное горение гомогенной водородно-воздушной смеси. Для описания кинетики химических превращений используется блок реакций окисления водорода из детального кинетического механизма окисления и горения нормальных углеводородов [17]. Значения коэффициентов полиномов для всех веществ взяты из работы [18]. Коэффициенты переноса $\lambda $ и ${{D}_{i}}$ рассчитывались по методике, описанной в [19].

Ниже данные, представленные в работе [20], используются для сравнения расчетных и измеренных значений скоростей турбулентного горения смесей разного состава в зависимости от интенсивности турбулентности. Максимальное физическое время расчета турбулентного пламени достигало 10^–3 с, т.е. было меньше характерного временнóго масштаба τ. Поэтому решения представляют собой индивидуальные реализации положения поверхности фронта турбулентного пламени. На поверхности фронта пламени наблюдаются (см. рис. 3) и мелкомасштабные, и крупномасштабные пространственные неоднородности, характеризуемые линейными размерами, которые и меньше, и больше толщины фронта пламени. Мелкомасштабная турбулентность изменяет (увеличивает) скорость обменных процессов (обмен массой и энергией) внутри фронта пламени, а крупномасштабная турбулентность изменяет форму пламени и увеличивает площадь его поверхности.

Определим истинную скорость турбулентного пламени u_t как скорость движения разделительной плоскости между свежей смесью и продуктами горения, т.е. как скорость условного фронта пламени. Несмотря на относительно короткое физическое время расчета (t < τ), истинная скорость турбулентного пламени почти постоянна (см. наклон кривой на рис. 5). Проверка показала, что одновременное удвоение числа расчетных точек по трем осям при сохранении размеров расчетной области (вместо 500 000 получается 4 000 000 ячеек) не изменяет значения истинной скорости u_t турбулентного пламени.

Если учесть тепловое расширение продуктов реакции, то по видимой скорости турбулентного пламени u_t можно определить истинную скорость его распространения U, используя соотношение

где m₀ и m₁ – число молей в свежей смеси и в продуктах реакции соответственно.

Отметим, что в редких случаях в расчетах получались сильно заниженные значения скорости (меньше нормальной скорости ламинарного пламени). Поскольку речь идет о синтетической турбулентности, вероятно, причина этого эффекта – необычные комбинации случайных чисел при нахождении функции y(t). При этом сама функция имела хотя и странную, но легко узнаваемую форму. Такие решения считались “паразитными” и не принимались во внимание.

На рис. 6 проведено, как и в [20], сравнение расчетных и измеренных зависимостей истинной скорости u_t распространения турбулентного пламени в смесях разного состава от пространственной среднеквадратичной пульсационной скорости $\bar {u}.$ Видно, что при увеличении турбулентности u_t возрастает. Из приведенного сравнения следует вывод об удовлетворительном качественном согласии результатов.

На рис. 7 проведено сравнение расчетных профилей температуры в ламинарном и турбулентном пламенах. По этим профилям можно оценить толщину пламени и сравнить ее с пространственными масштабами пульсаций скорости. Профиль температуры в турбулентном пламени более пологий, чем в ламинарном пламени, а само турбулентное пламя шире, чем ламинарное вследствие влияния мелкомасштабных обменных процессов.

На рис. 8 представлен пример расчетных профилей температуры и концентраций промежуточных и конечных продуктов реакции в турбулентном пламени. Несмотря на схожесть структур зон реакции в турбулентном и ламинарном пламенах, здесь все же есть некоторые количественные различия. Так, концентрации самых активных промежуточных продуктов реакции – гидроксила ОН, атомов Н и О – в турбулентном пламени оказываются ниже, чем в ламинарном. Это отчетливо видно на рис. 9, на котором представлены примеры расчетных зависимостей максимальных по зоне реакции концентраций гидроксила от среднеквадратичной пульсационной скорости $\bar {u}.$ Понижение концентрации активных центров является следствием увеличения интенсивности турбулентности, т.е. следствием ускорения обменных процессов. По косвенным данным такой вывод был сделан и в экспериментах, описанных в работе [21].

В заключение отметим еще один результат, который важен для понимания структуры зоны горения в турбулентном пламени: во всех проведенных расчетах поверхность фронта пламени была односвязной. Анализ расчетных распределений концентраций реагентов показал, что даже при максимальных интенсивностях турбулентности в поле течения не обнаруживаются “островки” свежей горючей смеси, окруженные продуктами горения. Это не согласуется с результатами расчетов в работе [10], где отмечалось наличие таких “островков”. По-видимому, здесь сказывается различие в методиках описания поля турбулентности или недостаточно длительное моделирование процесса распространения турбулентного пламени. Этот вопрос будет исследован дополнительно.

ЗАКЛЮЧЕНИЕ

Предложена методика трехмерного прямого численного моделирования распространения турбулентного пламени в газовых реагирующих смесях в поле синтетической стационарной, однородной и изотропной турбулентности. Методика основана на детальном кинетическом механизме горения многокомпонентной смеси, не содержит каких-либо подгоночных параметров и применена к расчету турбулентного горения водородно-воздушных смесей разного состава.

Сравнение результатов расчетов с экспериментальными данными показало, что между ними есть удовлетворительное качественное соответствие: согласно как расчету, так и эксперименту скорость турбулентного горения возрастает с увеличением интенсивности турбулентности. Расчеты показывают, что концентрации активных центров реакции – гидроксила ОН, атомов Н и О в турбулентном пламени меньше, чем в ламинарном, что также согласуется с экспериментом.

Дальнейшая работа будет направлена на применение предложенной методики к воздушным смесям углеводородных горючих, прежде всего метана. Важнейшая особенность турбулентного горения углеводородов, наблюдаемая экспериментально, – это погасание пламени при высоких интенсивностях турбулентности.

Работа выполнена при поддержке Российским научным фондом (проект 18-73-10196).