Химическая физика, 2019, T. 38, № 1, стр. 27-37
Прямое численное моделирование турбулентного горения водородно-воздушных смесей разного состава в двумерном приближенииВ. Я. Басевич, А. А. Беляев, С. М. Фролов, Ф. С. Фролов
В. Я. Басевич 1, А. А. Беляев 1, С. М. Фролов 1, 2, 3, *, Ф. С. Фролов 1, 3
1 Институт химической физики им. Н.Н. Семёнова Российской академии наук
Москва, Россия
2 Национальный исследовательский ядерный университет “МИФИ”
Москва, Россия
3 Федеральный научный центр Научно-исследовательский институт системных исследований
Российской академии наук
Москва, Россия
* E-mail: smfrol@chph.ras.ru
Поступила в редакцию 14.03.2018
После доработки 20.03.2018
Принята к публикации 14.03.2018
Аннотация
Предложена методика двумерного прямого численного моделирования распространения турбулентного пламени в газовых реагирующих смесях в условиях стационарной, однородной и изотропной турбулентности. Методика основана на детальном кинетическом механизме горения многокомпонентной смеси и не содержит каких-либо подгоночных параметров. Эта методика применена к расчету турбулентного горения водородно-воздушной смеси. Предложено условие, позволяющее сравнивать результаты двумерных расчетов – зависимости скорости распространения пламени от интенсивности турбулентности – с реальным “трехмерным” экспериментом. Полученное согласие расчетных и измеренных зависимостей подтверждает применимость предложенного условия. Рассмотрено влияние давления на скорость распространения пламени. Расчетные концентрации активных центров реакции – гидроксила ОН, атомов Н и О – в турбулентном пламени меньше, чем в ламинарном, что также согласуется с экспериментом.
ВВЕДЕНИЕ
Прямое численное моделирование турбулентного горения предполагает неэмпирический подход на основе применения детальных кинетических механизмов (ДКМ) и учета поля пульсационных скоростей. Первые работы по математическому моделированию турбулентного горения появились достаточно давно [1, 2], однако общепринятых методов математического описания физико-химических процессов в турбулентном пламени до сих пор в литературе нет [3–9]. Существующие подходы, как правило, включают всякого рода “замыкающие” гипотезы, основанные на экспериментальных наблюдениях. Поскольку область применимости той или иной гипотезы ограничена, такие полуэмпирические подходы неуниверсальны. Сегодня наиболее перспективным универсальным подходом к неэмпирическому теоретическому описанию турбулентного горения считают прямое численное моделирование (ПЧМ, в англоязычной литературе – Direct Numerical Simulation (DNS)), при котором в рассмотрение включены все основные особенности трехмерного турбулентного реагирующего течения с полным спектром турбулентных пульсаций скорости, с полным набором исходных, промежуточных и конечных химических компонентов с их индивидуальными термохимическими свойствами и свойствами молекулярного переноса, а также с адекватными граничными условиями [7–9]. Численное интегрирование определяющих уравнений течения проводят на расчетных сетках, обеспечивающих пространственное разрешение турбулентных вихрей колмогоровского масштаба, используя схемы высокого порядка аппроксимации.
Несмотря на значительные успехи в развитии ПЧМ турбулентного горения, на этом пути есть еще много проблем. Так, при сравнении самых “свежих” опубликованных решений для скорости распространения турбулентного пламени, полученных с помощью ПЧМ (см., например, [10]), с известными экспериментальными данными выясняется, что их точность все еще заметно уступает точности, достигаемой при решении задач ламинарного горения. В теоретических работах в качестве причин расхождения результатов расчетов и экспериментов часто называют неучтенные эффекты, вызванные неопределенными граничными условиями в реальных экспериментах. Кроме того, применение ПЧМ к реальным турбулентным пламенам требует очень больших вычислительных ресурсов, в связи с чем в расчетах все еще приходится использовать различные упрощения, влияющие на точность решения.
В данной работе используется альтернативный подход к ПЧМ турбулентного горения гомогенной газовой смеси, который ранее был предложен в работе [11]. Вместо численного решения всех уравнений, определяющих распространение турбулентного пламени в реагирующем газе, в [11] было предложено решать только уравнения переноса скалярных величин – концентраций реагентов и энергии в искусственном (“синтетическом”) поле турбулентности, характеризуемом заданной (постоянной) среднеквадратичной интенсивностью пульсаций скорости и заданными (постоянными) интегральными пространственными и временны́ми масштабами. При этом считается, что распространение пламени не влияет на характеристики синтетического поля турбулентности в предпламенной зоне.
В отличие от [11] в данной работе химический процесс описывается не одной глобальной реакцией, а детальным кинетическим механизмом с участием активных центров. Первые результаты по такому моделированию турбулентного горения представлены нами в работе [12] для одного конкретного состава водородно-воздушной смеси в простейших условиях стационарной, однородной и изотропной турбулентности. Цель данной работы – совершенствование методики из работы [12] и ее применение для расчетов распространения турбулентного пламени в водородно-воздушных гомогенных смесях разного состава и при разных давлениях, а также сравнение результатов расчетов с известными экспериментальными данными.
МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ
Система уравнений, определяющих распространение турбулентного пламени, основана на уравнениях Навье–Стокса в трехмерной постановке, уравнении сохранения энергии и уравнениях неразрывности для всех химических компонентов в смеси идеальных газов [13]:
(1)
$\rho \frac{{\partial {{Y}_{i}}}}{{\partial t}} + \rho {\mathbf{v}}\nabla {{Y}_{i}} = {{\omega }_{i}} - \nabla (\rho {{Y}_{i}}{{{\mathbf{V}}}_{i}}),\;\;i = 1,...,N,$где $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) сводится к виду
(2)
$\rho \frac{{\partial e}}{{\partial t}} + \rho {\mathbf{v}}\nabla e = - \nabla (\lambda \nabla T),$где $\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{'}}){\text{:}}$
Здесь и далее прописные буквы означают среднее значение, а штрихи относятся к пульсациям скорости. Тогда (3) можно преобразовать к виду
(4)
$\nabla (\rho {\mathbf{v}}) = \nabla (\rho {\mathbf{V}} + \rho {\mathbf{v}}{\text{'}}) = \nabla (\rho {\mathbf{V}}) + \nabla (\rho {\mathbf{v}}{\text{'}}) = 0.$Основываясь на первом предположении об однородной и изотропной турбулентности, предположим далее, что
и, следовательно,
В области течения с простейшей геометрией (см. допущение (1)) всегда существует усредненное направление распространения турбулентного пламени, и можно устремить одну ось пространственных координат в этом направлении. Тогда (6) можно проинтегрировать и получить
где ${\mathbf{B}}$ – постоянный поток вещества. С учетом (7) систему (2) можно привести к виду
(8)
$\begin{gathered} \rho \frac{{\partial e}}{{\partial t}} + {\mathbf{B}}\nabla e = - \nabla (\lambda \nabla T) - \nabla (\rho {\mathbf{v}}{\text{'}}u), \\ \rho \frac{{\partial {{Y}_{i}}}}{{\partial t}} + {\mathbf{B}}\nabla {{Y}_{i}} = {{\omega }_{i}} - \nabla ({{D}_{i}}\rho \nabla {{Y}_{i}}) - \nabla (\rho {\mathbf{v}}{\text{'}}{{Y}_{i}}), \\ i = 1,...,N. \\ \end{gathered} $Последние члены в обоих уравнениях системы (8) определяют перенос вещества и энергии с турбулентными пульсациями скорости. Если вектор ${\mathbf{B}}$ известен, то производные по времени в обоих уравнениях системы (8) можно положить равными нулю, чтобы получить структуру стационарного турбулентного пламени. Однако, как правило, значение вектора ${\mathbf{B}}$ заранее неизвестно. Поэтому, чтобы получить его значение, необходимо решить нестационарную систему (8), описывающую движение волны горения до достижения некоторой постоянной скорости. В этом случае конвективные члены в системе (8) можно опустить, а система уравнений примет вид
(9)
$\begin{gathered} \rho \frac{{\partial e}}{{\partial t}} = - \nabla (\lambda \nabla T) - \nabla (\rho {\mathbf{v}}{\text{'}}e), \\ \rho \frac{{\partial {{Y}_{i}}}}{{\partial t}} = {{\omega }_{i}} - \nabla ({{D}_{i}}\rho \nabla {{Y}_{i}}) - \nabla (\rho {\mathbf{v}}{\text{'}}{{Y}_{i}}), \\ i = 1,...,N. \\ \end{gathered} $Наконец, можно упростить и последние члены в двух уравнениях (9), используя (5) и предполагая, что
(10)
${{\nabla }_{x}}(\rho u{\text{'}}) = {{\nabla }_{y}}(\rho v{\text{'}}) = {{\nabla }_{z}}(\rho w{\text{'}}) = 0,$где ${{\nabla }_{x}} = \frac{\partial }{{\partial x}},$ ${{\nabla }_{y}} = \frac{\partial }{{\partial y}}$ и ${{\nabla }_{z}} = \frac{\partial }{{\partial z}}.$ Интегрирование уравнений (10) дает
(11)
$\begin{gathered} \rho u{\text{'}} = {{\rho }_{0}}u_{0}^{'}, \\ \rho v{\text{'}} = {{\rho }_{0}}v_{0}^{'}, \\ \rho w{\text{'}} = {{\rho }_{0}}w_{0}^{'}. \\ \end{gathered} $После подстановки соотношений (11) в систему (9) окончательно получим
(12)
$\begin{gathered} \rho \frac{{\partial e}}{{\partial t}} = - \nabla (\lambda \nabla T) - {{\rho }_{0}}{{c}_{v}}{\mathbf{v}}_{0}^{'}\nabla T, \\ \rho \frac{{\partial {{Y}_{i}}}}{{\partial t}} = {{\omega }_{i}} - \nabla ({{D}_{i}}\rho \nabla {{Y}_{i}}) - {{\rho }_{0}}{\mathbf{v}}_{0}^{'}\nabla {{Y}_{i}}, \\ i = 1,...,N. \\ \end{gathered} $Таким образом, как и в [11, 12], вместо численного решения системы (1) предлагается решать систему (12), состоящую только из уравнений переноса скалярных величин – концентраций всех N реагентов и энергии на заданном синтетическом поле турбулентности, характеризуемом среднеквадратичной интенсивностью пульсаций скорости v' и заданными интегральными пространственными и временны́ми масштабами. В качестве первого шага можно ограничиться решением двумерной задачи, хотя в двумерном приближении поверхность пламени меньше, чем в трехмерном, а двумерная турбулентность воздействует на распространение турбулентного пламени слабее, чем трехмерная. В двумерном приближении в системе координат (x, y) уравнения (12) с учетом калорического уравнения состояния примут вид
(13)
$\begin{gathered} \rho {{c}_{p}}\frac{{\partial T}}{{\partial t}} = \sum\limits_{i = 1}^N {{{h}_{i}}{{\omega }_{i}}} + \left( {\frac{\partial }{{\partial x}}\lambda \frac{{\partial T}}{{\partial x}} + \frac{\partial }{{\partial y}}\lambda \frac{{\partial T}}{{\partial v}}} \right) - \\ - \,\,{{\rho }_{0}}{{c}_{p}}\left( {u_{0}^{'}\frac{{\partial T}}{{\partial x}} + v_{0}^{'}\frac{{\partial T}}{{\partial y}}} \right), \\ \rho \frac{{\partial {{Y}_{i}}}}{{\partial t}} = {{\omega }_{i}} + \left( {\frac{\partial }{{\partial x}}{{D}_{i}}\rho \frac{{\partial {{Y}_{i}}}}{{\partial x}} + \frac{\partial }{{\partial y}}{{D}_{i}}\rho \frac{{\partial {{Y}_{i}}}}{{\partial y}}} \right) - \\ - \,\,{{\rho }_{0}}\left( {u_{0}^{'}\frac{{\partial {{Y}_{i}}}}{{\partial x}} + v_{0}^{'}\frac{{\partial {{Y}_{i}}}}{{\partial y}}} \right),\,\,\,\,i = 1,...,N, \\ {{p}_{0}} = \rho {{R}^{0}}T\sum\limits_{i = 1}^N {\frac{{{{Y}_{i}}}}{{{{\mu }_{i}}}}} . \\ \end{gathered} $Система (13) состоит из N + 2 уравнений для N + 2 переменных ($N$ веществ с массовой долей ${{Y}_{i}}$, температурой $T$ и плотностью $\rho $). Для замыкания системы используем следующие соотношения для ${{c}_{p}},$ $\lambda $ и ${{D}_{i}}{\text{:}}$
(14)
$\begin{gathered} {{c}_{p}} = \sum\limits_{i = 1}^N {{{c}_{{p,i}}}{{Y}_{i}}} , \\ {{c}_{{p,i}}} = \frac{{{{R}^{0}}}}{{{{\mu }_{i}}}}\left( {{{a}_{1}} + {{a}_{2}}T + {{a}_{3}}{{T}^{2}} + {{a}_{4}}{{T}^{3}} + {{a}_{5}}{{T}^{4}}} \right), \\ \lambda = \lambda (T,{{Y}_{1}},...,{{Y}_{N}},{{\mu }_{1}},...,{{\mu }_{N}},{{c}_{{p,1}}},...,{{c}_{{p,N}}}), \\ {{D}_{i}} = {{D}_{i}}(T,p,{{\mu }_{i}},{{\mu }_{{in}}}), \\ \end{gathered} $где a1, a2, a3, a4 и a5 – коэффициенты полиномов, а индекс “in” относится к инертному разбавителю (азот).
Для скорости химических реакций используется соотношение [13]
где $\nu _{{i,k}}^{ + }$ и $\nu _{{i,k}}^{ - }$ – стехиометрические коэффициенты для i-тых веществ, являющихся реагентами и продуктами k-й реакции, соответственно; ${{A}_{k}}$ – предэкспоненциальный множитель k-й реакции; ${{\alpha }_{k}}$ – показатель степени, определяющий температурную зависимость предэкспоненциального множителя k-й реакции; ${{E}_{k}}$ – энергия активации k-й реакции; $M$ – общее количество химических реакций и ${{X}_{j}}$ – мольная доля j-го вещества.
Уравнения (13) содержат параметры стационарной, однородной и изотропной турбулентности реагирующего газа (две компоненты вектора ${\mathbf{v}}_{0}^{'} = (u_{0}^{'},v_{0}^{'})$). Компоненты вектора пульсационной скорости можно получить статистическим разыгрыванием по методу Монте-Карло, предположив, что мгновенные местные компоненты вектора пульсационной скорости удовлетворяют нормальному распределению Гаусса φ, а вихревая структура турбулентности описывается экспоненциально затухающими пространственными (R' и R") и временной (R) корреляционными функциями:
(15)
$\varphi (u{\text{'}}) = \frac{1}{{\sqrt {2\pi } \sigma }}\exp \left[ { - \frac{{{{{(u{\text{'}} - \bar {u})}}^{2}}}}{{2{{\sigma }^{2}}}}} \right],$где $u{\text{'}}$ – мгновенная составляющая вектора пульсационной скорости, $\bar {u}$ – длина вектора средней пульсационной скорости, $\sigma $ – среднеквадратичное отклонение пульсаций скорости, rk – расстояние в физическом пространстве, L' и L" – масштабы турбулентности в направлениях y и x, соответственно, τ – временнóй масштаб Лагранжа. Несмотря на допущение об изотропности поля турбулентности (допущение 1, см. выше) здесь для общности введены два пространственных масштаба: L' и L". В соответствии с допущением о стационарности и однородности поля турбулентности масштабы τ, L' и L" приняты постоянными. Граничные и начальные условия для системы (13) рассмотрены ниже для конкретной расчетной области.
МЕТОДИКА РАСЧЕТОВ
На рис. 1 показана схема простейшей расчетной области с выбранной системой координат (x, y). Расчетная область представляет собой прямоугольник, в котором левая и правая границы – непроницаемые стенки со скольжением потока, а нижняя и верхняя границы удалены от пламени на достаточно большое расстояние, чтобы обеспечить постоянство давления в системе за все время расчета. Начальное положение плоского фронта пламени в момент времени t = t0= 0 показано нижней горизонтальной штриховой линией. Пламя распространяется снизу–вверх: выше штриховой линии находится свежая горючая смесь, а ниже – продукты горения. Здесь и далее при t > 0 положение фронта пламени определяется как геометрическое место точек, в которых температура равна среднему арифметическому между начальной температурой свежей смеси T0 и термодинамически равновесной температурой продуктов сгорания T1: Тm = (T0+ T1)/2.
Мгновенное положение фронта пламени в момент времени t = ti > 0 показано сплошной линией, а усредненное положение фронта пламени – верхней горизонтальной штриховой линией, которая проведена таким образом, чтобы площади областей, отсекаемых сплошной кривой выше и ниже штриховой линии, были одинаковы.
Задача о распространении турбулентного пламени решалась с помощью численного интегрирования системы (13) с дополнительными соотношениями методом переменных направлений [14]. Для решения уравнений по каждому из направлений (x, y) применялась неявная разностная схема с равномерным шагом по пространству и линеаризацией нелинейных “источниковых” слагаемых на верхнем слое. Эта схема имеет первый порядок точности по времени и пространству. Расчетная область имеет размер 0.5 × 4 см2. Все расчетные ячейки имеют форму квадрата размером 0.005 × 0.005 см2. Полное количество расчетных ячеек равно 80 000. Шаг интегрирования по времени изменяется в зависимости от числа итераций, но не превышает 10–6 с.
Важная операция – моделирование синтетического поля турбулентности ${\mathbf{v}}_{0}^{'} = (u_{0}^{'},v_{0}^{'})$ в (13). Моделирование полей пульсационных скоростей $u{\text{'}}$ = $u_{x}^{'}$(x, y, t) и $v{\text{'}}$ = $u_{y}^{'}$(x, y, t) проводится путем разыгрывания их возможных значений при заданных характеристиках турбулентности $\bar {u},$ L', L" и τ. Соответствующая методика описана в Приложении. Подчеркнем, что принятая методика получения синтетического поля турбулентности не требует задания ни спектра турбулентности, ни колмогоровского масштаба. Мгновенная длина вектора пульсационной скорости, u, определяется по значениям $u_{x}^{'}$ и $u_{y}^{'}{\text{:}}$
Ниже представлены используемые здесь начальные и граничные условия, соответственно:
(17)
$\begin{gathered} t = {{t}_{0}} = 0,\,\,\,\,y > y(t = {{t}_{0}}){\text{:}}\,\,{{Y}_{j}} = {{Y}_{{j0}}},T = {{T}_{0}}, \\ j = 1,...,N, \\ y < y(t = {{t}_{0}}){\text{:}}\,\,{{Y}_{j}} = {{Y}_{{j{\text{1}}}}},\,\,\,\,T = {{T}_{{\text{1}}}},j,...,N; \\ \end{gathered} $(18)
$\begin{gathered} х = 0{\text{:}}\,\,\frac{{\partial T}}{{\partial x}} = 0,\,\,\,\,\frac{{\partial {{Y}_{j}}}}{{\partial x}} = 0,\,\,\,\,j = {\text{1}},...,N, \\ х = {{L}_{x}}{\text{:}}\,\,\frac{{\partial T}}{{\partial x}} = 0,\,\,\,\,\frac{{\partial {{Y}_{j}}}}{{\partial x}} = 0,\,\,\,\,j = {\text{1}},...,N, \\ y = 0{\text{:}}\,\,\quad\frac{{\partial T}}{{\partial y}} = 0,\,\,\,\,\frac{{\partial {{Y}_{j}}}}{{\partial y}} = 0,\,\,\,\,j = {\text{1}},...,N, \\ y = {{L}_{y}}{\text{:}}\,\,\quad\frac{{\partial T}}{{\partial y}} = 0,\,\,\,\,\frac{{\partial {{Y}_{j}}}}{{\partial y}} = 0,\,\,\,\,j = {\text{1}},...,N. \\ \end{gathered} $РЕЗУЛЬТАТЫ РАСЧЕТОВ
Как и в [12], в данной работе рассматривается турбулентное горение гомогенной водородно-воздушной смеси. Для описания кинетики химических превращений используется блок реакций окисления водорода из ДКМ окисления и горения нормальных углеводородов [15]. Значения коэффициентов полиномов для всех веществ взяты из [16]. Коэффициенты переноса λ и ${{D}_{i}}$ рассчитывались по методике, описанной в [17].
Прежде чем решать задачу о распространении турбулентного пламени мы провели численное интегрирование системы (13) с начальными и граничными условиями (17), (18) без синтетической турбулентности, т.е. с нулевыми пульсационными скоростями ${\mathbf{v}}_{0}^{'} = (u_{0}^{'},v_{0}^{'}) = 0,$ чтобы проверить применимость методики к описанию распространения плоского ламинарного пламени со скоростью un. Расчеты проведены для условий с Т0 = 293 К и p = 1 ата. Объемное содержание водорода в смеси изменялось от 9.09% до 23.09%. На рис. 2 результаты этих расчетов представлены символами 1. Полученные результаты сравниваются на рис. 2 c экспериментальными данными [3, 18–23], а также с расчетами по программе [24] для более широкого интервала объемного содержания водорода: от 6% до 75%. Отметим, что в работе [18] представлены экспериментальные данные как по ламинарным, так и по турбулентным пламенам, полученные для смесей с объемным содержанием Н2 от 9.09% до 23.09%. Ниже именно эти данные [18] используются для сравнения расчетных и измеренных скоростей турбулентного горения смесей разного состава от интенсивности турбулентности. Для водородно-воздушных смесей с содержанием водорода Н2 < 9.0% получить cтационарное решение для плоского ламинарного пламени не удалось.
На рис. 3 показаны примеры мгновенной формы фронта турбулентного пламени в различные моменты времени после начала расчета. Максимальное физическое время расчета турбулентного пламени достигало 10–3 с, т.е. было меньше характерного временнóго масштаба τ. Кривые на рис. 3 представляют собой индивидуальные реализации положения фронта турбулентного пламени. На фронте пламени наблюдаются и мелкие, и крупные пространственные неоднородности. Как будет показано ниже, здесь присутствуют линейные размеры, которые и меньше, и больше толщины фронта пламени. Мелкомасштабная турбулентность изменяет (увеличивает) скорость обменных процессов (обмен массой и энергией) внутри фронта пламени, а крупномасштабная турбулентность изменяет форму пламени и увеличивает площадь его поверхности.
Несмотря на относительно короткое физическое время расчета (t < τ), полученное значение скорости распространения усредненного турбулентного пламени почти постоянно во времени (рис. 4). По положениям фронта пламени, усредненным по нескольким реализациям (обычно 5–7 реализаций), можно определить видимую скорость его распространения U при данной пульсационной скорости. Если учесть тепловое расширение продуктов реакции, можно определить скорость распространения турбулентного пламени:
где m0 и m1 – число молей в свежей смеси и в продуктах реакции соответственно.
На рис. 5 сравниваются расчетные и измеренные зависимости скорости распространения турбулентного пламени ut в смесях разного состава от среднеквадратичной пульсационной скорости $\bar {u}.$ Отметим одно важное обстоятельство. На рис. 5, как и у авторов экспериментальной работы [18], по оси абсцисс отложена величина $\bar {u} = \sqrt 3 u{\text{'}}$ – среднеквадратичное значение пульсационной скорости, где u' – проекция $\bar {u}$ на ось y. Согласно (4), для корректного сравнения двумерных расчетов с реальным “трехмерным” экспериментом вдоль оси абсцисс должна быть отложена величина $\bar {u} = \sqrt 2 u{\text{'}},$ т.е. для расчетных значений пульсационной скорости линейный масштаб должен быть уменьшен в отношении $\sqrt 2 :\sqrt 3 \approx 0.815.$ Именно это условие позволяет сравнивать результаты двумерных расчетов с “трехмерным” экспериментом. Из приведенного сравнения расчетных и измеренных результатов следует вывод об их удовлетворительном согласии.
На рис. 6 сравниваются расчетные профили температуры в ламинарном и турбулентном пламенах (мгновенная реализация). По этим профилям можно оценить толщину пламени и сравнить ее с пространственными масштабами пульсаций скорости. Профиль температуры в турбулентном пламени более пологий, чем в ламинарном пламени, а само турбулентное пламя шире, чем ламинарное, вследствие влияния мелкомасштабных обменных процессов.
На рис. 7 представлен пример расчетных профилей температуры и концентраций промежуточных и конечных продуктов реакции в турбулентном пламени (мгновенная реализация). Несмотря на схожесть структуры зон реакции турбулентного и ламинарного пламен, здесь все же есть некоторые количественные различия. Так, концентрации самых активных промежуточных продуктов реакции – гидроксила ОН, атомов Н и О в турбулентном пламени оказываются ниже, чем в ламинарном. Это отчетливо видно из рис. 8, на котором представлены примеры расчетных зависимостей максимальных по зоне реакции концентраций гидроксила от среднеквадратичной пульсационной скорости $\bar {u}.$ Понижение концентрации активных центров есть следствие увеличения интенсивности турбулентности, т.е. ускорения обменных процессов. По косвенным данным такой вывод был сделан и в экспериментах [25].
Хорошо известно, что с повышением давления скорость распространения ламинарного пламени уменьшается, а скорость распространения турбулентного пламени либо остается постоянной, либо несколько увеличивается [5]. Последнее обнаруживается и в наших расчетах распространения турбулентного пламени при разных давлениях. На рис. 9 сравниваются расчетные зависимости для скоростей распространения ламинарного и турбулентного пламен.
В заключение отметим один результат, который важен для понимания структуры зоны горения в турбулентном пламени: во всех проведенных расчетах поверхность фронта пламени была односвязной. Даже при максимальных интенсивностях турбулентности расчеты не показали наличия “островков” свежей горючей смеси, окруженных продуктами горения, как это было описано в [10]. По-видимому, здесь сказываются различия в методиках описания синтетического поля турбулентности или недостаточно длительное моделирование процесса распространения турбулентного пламени. Этот вопрос будет исследоваться дополнительно.
ЗАКЛЮЧЕНИЕ
Предложена методика двумерного ПЧМ распространения турбулентного пламени в газовых реагирующих смесях в условиях стационарной, однородной и изотропной синтетической турбулентности. Методика основана на ДКМ горения многокомпонентной смеси и не содержит каких-либо подгоночных параметров. Она применена к расчету турбулентного горения водородно-воздушных смесей разного состава при разных начальных давлениях. Предложено условие, позволяющее сравнивать результаты двумерных расчетов с реальным экспериментом.
Сравнение результатов расчетов с экспериментальными данными показало, что между ними есть удовлетворительное качественное согласие: и в расчете, и в эксперименте скорость турбулентного горения возрастает с увеличением интенсивности турбулентности. Правильно отражается и зависимость скорости распространения турбулентного пламени от давления. Кроме того, расчеты показывают, что концентрации активных центров реакции – гидроксила ОН, атомов Н и О в турбулентном пламени меньше, чем в ламинарном, что также согласуется с экспериментом. Полученное качественное согласие подтверждает применимость предложенного условия для сравнения двумерных расчетов с реальным экспериментом. Дальнейшая работа будет направлена на расширение предложенной методики на трехмерный случай и учет вихревой структуры турбулентности.
Работа выполнена при частичной поддержке Российским фондом фундаментальных исследований (проект офи-м № 16-29-01065), а также за счет субсидии, выделенной ИХФ РАН на выполнение государственного задания ФАНО России (темы 0082-2016-0011, AAA-A17-117040610346-5 и 0082-2014-0004, AAA-A17-117040610283-3).
Приложение
Моделирование поля синтетической турбулентности
Геометрические характеристики поля: х – ширина (масштаб пульсаций L"), у – направление движения пламени (масштаб пульсаций L'), Δxk и Δt – пространственные и временны́е шаги интегрирования, k = 1, 2 (x1 = x, x2= y).
1. Вычисление ux(х = 0,…, N, у = 0, t = 0):
– в точке (х = 0, у = 0, t = 0) в соответствии с заданным значением $\bar {u}$ генератор случайных чисел выдает ux(0, 0, 0); здесь и далее перед каждым вызовом генератора случайных чисел необходимо проверить среднее значение (СЗ): если |СЗ| > u, заданной в начальных условиях, то следует положить |СЗ| = u;
– в точке (х = Δx1, 0, 0) вычисляются значение R", среднее значение пульсационной скорости ux (0, 0, 0)R", дисперсия ${{\bar {u}}^{2}}\left( {{\text{1}}--{{R}^{{''2}}}} \right)$ и разыгрывается значение пульсационной скорости ux(Δx1, 0, 0);
– аналогично через интервалы Δxi до конца расчетного отрезка х определяются ux(1, 0, 0), ux(2, 0, 0), …, ux(N, 0, 0).
2. Вычисление uy(х = 0, …, N, у = 0, t = 0):
– в точке (х = 0, у = 0, t = 0) в соответствии с заданным $\bar {u}$ генератор случайных чисел выдает uу(0, 0, 0);
– в точке (х = Δx1, 0, 0) вычисляются значение R", среднее значение пульсационной скорости uу(0, 0, 0)R", дисперсия ${{\bar {u}}^{2}}\left( {{\text{1}}--{{R}^{{''2}}}} \right)$ и разыгрывается значение пульсационной скорости uу(Δx1, 0, 0);
– аналогично через интервалы Δxi до конца расчетного отрезка х определяются uу(1, 0, 0), uу(2, 0, 0), …, uу(N, 0, 0).
3. Вычисление uy(х = 0, …, N, у = 1, t = 0):
– по полученным в п.1 значениям ux в точках 1, 2, …, N вычисляются производные ${{\partial {{u}_{x}}} \mathord{\left/ {\vphantom {{\partial {{u}_{x}}} {\partial x}}} \right. \kern-0em} {\partial x}}$ по всей оси х(0 – N, 0, 0);
– по уравнению неразрывности определяются значения производных ${{\partial {{u}_{у }}} \mathord{\left/ {\vphantom {{\partial {{u}_{у }}} {\partial y}}} \right. \kern-0em} {\partial y}}$ по всей оси х (0 – N, 0, 0);
– в точке (х = 0, 1, 0) вычисляется
– в точке (х = 0, у = 1, t = 0) вычисляются значение R', среднее значение пульсационной скорости [〈uу(0, 1, 0)〉]R', дисперсия ${{\bar {u}}^{2}}\left( {{\text{1}}--{{R}^{{'2}}}} \right)$ и разыгрывается вероятное значение пульсационной скорости uу(0, 1, 0);
– аналогично через интервалы Δxi до конца расчетного отрезка х определяются uу(1, 1, 0), uу(2, 1, 0), …, uу(N, 1, 0);
– для всех uу(N, 1, 0) вычисляется получившееся среднеквадратичное значение 〈uу(0, …, N, 1, 0)〉, определяется корректирующий коэффициент
находятся скорректированные значения uу(0, …, N, 1, 0) = uу(0, …, N, 1, 0)/au.
4. Вычисление ux(х = 0, …, N, у = 1, t = 0):
– по полученным в пп. 1, 3 значениям uy в точках 1, 2, …, N вычисляются производные ∂uy/∂у по всей оси х(0, …, N, 0, 0);
– по уравнению неразрывности определяются значения производных ∂uх/∂х по всей оси х(0, …, N, 1, 0);
– в точке (х = 0, у = 1, t = 0) в соответствии с заданным $\bar {u}$ генератор случайных чисел выдает ux(0, 1, 0);
– в точке (х = 1, 1, 0) вычисляется
– в точке (х = 1, у = 1, t = 0) вычисляется значение R", среднее значение пульсационной скорости [〈uх(1, 1, 0)〉]R", дисперсия ${{\bar {u}}^{2}}\left( {{\text{1}}--{{R}^{{''2}}}} \right)$ и разыгрывается вероятное значение пульсационной скорости uх(1, 1, 0);
– аналогично через интервалы Δxi до конца расчетного отрезка х определяются uх(1, 1, 0), uх(2, 1, 0), …, uх(N, 1, 0);
– для всех uх(N, 1, 0) вычисляется получившееся среднеквадратичное значение 〈uх(0, …, N, 1, 0)〉, определяется корректирующий коэффициент
находятся скорректированные значения uх(0, …, N, 1, 0) = uх(0, …, N, 1, 0)/aux.
5. Вычисление ux(х = 0, …, N, у = 1, …, N, t = 0):
– как в пп.1–4 .
6. Вычисление ux(х = 0, ..., N, у = 1, …, N, t = 0, …, tend) и uу(х = 0, …, N, у = 1, …, N, t = 0, …, tend):
– как изменяются значения uх и uу при переходе от оси х (0, …, N, 0, 0) к оси х (0, …, N, 1, 0), так они изменяются и при переходе от оси х (0, …, N, 0, 0) к оси х (0, …, N, 0, 1), но с использованием коэффициента корреляции τ. Эта методика сохраняется на протяжении всего расчета: от t = 0 до tend.
Проверка алгоритма. Получаемые поля пульсационных скоростей следует проверить на соблюдение задаваемых характеристик турбулентности $\bar {u},$ соответствие распределения пульсационных скоростей нормальному закону Гаусса, соблюдение масштабов L', L" и τ.
Примечание. Из предположения, что для турбулентных пульсаций выполняются условия (10), следует связь “холодных” (индекс “0”) и “горячих” (Т > Т0) пульсаций скорости:
где uk, 0 – проекции мгновенной пульсационной скорости на ось хk при начальной плотности. Этим уравнением определяется изменение интенсивности турбулентности в процессе горения, а разыгрывание пульсаций скорости требуется только при начальной температуре.
Согласно алгоритму моделирования синтетического поля турбулентности значение среднеквадратичной пульсационной скорости $\bar {u}$ постоянно во времени, а для всего расчетного пространства постоянны характеристики турбулентности L', L" и τ.
Список литературы
Damkoler G. // Z. Elektrochem. 1940. B. 46. S. 601.
Щелкин К.И. // ЖТФ. 1943. Т. 13. С. 520.
Lewis B., Elbe G. Combustion, Flames and Explosions of Gases. Orlando: Acad. Press. 1987.
Соколик А.С. Самовоспламенение, пламя и детонация в газах. М.: Изд-во АН СССР, 1960.
Щетинков Е.С. Физика горения газов. М.: Наука, 1965.
Bernard J.A., Bradley J.N. Flame and Combustion. L.–N.Y.: Chapman. Hall, 1985.
Echekki T., Chen J.H. // Combust and Flame. 2003. V. 134. P. 169.
Bell J.B., Day M.S., Grcar J.F. // Proc. Combust. Inst. 2002. V. 29. P. 1987.
Bell J.B., Cheng R.K., Day M.S., Shepherd I.G. // Ibid. 2006. V. 31. P. 1309.
Aspden A.J., Day M.S., Bell J.B. // Combust and Flame. 2016. V. 166. P. 266.
Басевич В.Я., Володин В.П., Когарко С.М., Перегудов Н.И. // Хим. физика. 1982. Т. 1. № 8. С. 1130.
Басевич В.Я., Беляев А.А., Фролов С.М., Басара Б. // Горение и взрыв. 2017. Т. 10. № 1. С. 4.
Вильямс Ф.А. Теория горения. М.: Наука, 1971.
Годунов С.К., Рябенький В.С. Разностные схемы. М.: Наука, 1977.
Басевич В.Я., Беляев А.А., Посвянский В.С., Фролов С.М. // Хим. физика. 2013. Т. 32. № 4. С. 87.
Thermodynamic Data at the Web Site of the Laboratory for Chemical Kinetics. Alexander Burcat’s. Ideal Gas Thermodynamic Data in Polynomial Form for Combustion and Air Pollution Use; http://garfield.chem.elte.hu/ Burcat/burcat.html.
Рид Р., Праусниц Дж., Шервуд Т. Свойства газов и жидкостей. Л.: Химия, 1982.
Карпов В.П., Северин Е.С. // Физика горения и взрыва. 1980. Т. 16. № 1. С. 45.
Козаченко Л.С. Дис. … д-ра физ.-мат. наук. М.: ИХФ АН СССР, 1954.
Manton J., Milliken B.B. // Proc. Gas Dynamics Symp. (Aerothermochem.) Northwestern Univ., 1956. P. 151.
Andrews G.E., Bradley D. // Combust and Flame. 1973. V. 20. P. 77.
Iijima T., Takeno T. // Ibid. 1986. V. 65. P. 35.
Dowdy D.R., Smith D.B., Taylor S.C., Williams A. // Proc. Combust. Inst. 1990. V. 23. P. 325.
Беляев А.А., Посвянский В.С. // Алгоритмы и программы. Информ. бюлл. гос. фонда алгоритмов и программ СССР. 1985. № 3. С. 35.
Басевич В.Я., Когарко С.М. // Физика горения и взрыва. 1985. Т. 21. № 5. С. 12.
Дополнительные материалы отсутствуют.
Инструменты
Химическая физика