Химическая физика, 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

Полный текст (PDF)

Аннотация

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

Ключевые слова: прямое численное моделирование, турбулентное горение, детальный кинетический механизм, водород.

ВВЕДЕНИЕ

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

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

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

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

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

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

$\frac{{\partial \rho }}{{\partial t}} + \nabla (\rho {\mathbf{v}}) = 0,$
$\rho \frac{{\partial {\mathbf{v}}}}{{\partial t}} + \rho {\mathbf{v}}\nabla {\mathbf{v}} = - \nabla p + \rho \sum\limits_{i = 1}^N {{{Y}_{i}}{{{\mathbf{f}}}_{i}}} ,$
$\rho \frac{{\partial e}}{{\partial t}} + \rho {\mathbf{v}}\nabla e = - \nabla {\mathbf{q}} - {\mathbf{P}}{\text{:}}\,\,(\nabla {\mathbf{v}}) + \rho \sum\limits_{i = 1}^N {{{Y}_{i}}{{{\mathbf{f}}}_{i}}} {{{\mathbf{V}}}_{i}},$
(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,$
$p = \rho {{R}^{0}}T\sum\limits_{i = 1}^N {\frac{{{{Y}_{i}}}}{{{{\mu }_{i}}}}} ,$
$e = \sum\limits_{i = 1}^N {{{h}_{i}}{{Y}_{i}} - \frac{p}{\rho },} $
${{h}_{i}} = h_{i}^{0} + \int\limits_{{{T}_{0}}}^T {{{c}_{{p,i}}}dT} ,\;\;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) сводится к виду

$\frac{{\partial \rho }}{{\partial t}} + \nabla (\rho {\mathbf{v}}) = 0,$
(2)
$\rho \frac{{\partial e}}{{\partial t}} + \rho {\mathbf{v}}\nabla e = - \nabla (\lambda \nabla T),$
$\rho \frac{{\partial {{Y}_{i}}}}{{\partial t}} + \rho {\mathbf{v}}\nabla {{Y}_{i}} = {{\omega }_{i}} - \nabla ({{D}_{i}}\rho \nabla {{Y}_{i}}),\,\,i = 1,...,N,$

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

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

(3)
$\nabla (\rho {\mathbf{v}}) = 0.$

Поскольку нас интересует именно такое стационарное решение, уравнение (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{:}}$

${\mathbf{v}} = {\mathbf{V}} + {\mathbf{v}}'.$

Здесь и далее прописные буквы означают среднее значение, а штрихи относятся к пульсациям скорости. Тогда (3) можно преобразовать к виду

(4)
$\nabla (\rho {\mathbf{v}}) = \nabla (\rho {\mathbf{V}} + \rho {\mathbf{v}}{\text{'}}) = \nabla (\rho {\mathbf{V}}) + \nabla (\rho {\mathbf{v}}{\text{'}}) = 0.$

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

(5)
$\nabla (\rho {\mathbf{v}}{\text{'}}) = 0,$

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

(6)
$\nabla (\rho {\mathbf{V}}) = 0.$

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

(7)
$\rho {\mathbf{V}} = {{\rho }_{0}}{{{\mathbf{V}}}_{0}} = {\mathbf{B}},$

где ${\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]

$\begin{gathered} {{\omega }_{i}} = {{\mu }_{i}}\sum\limits_{k = 1}^M {(\nu _{{i,k}}^{ - } - \nu _{{i,k}}^{ + }){{A}_{k}}{{T}^{{{{\alpha }_{k}}}}}} \times \\ \times \,\,\exp \left\{ { - ({{{{E}_{k}}} \mathord{\left/ {\vphantom {{{{E}_{k}}} {{{R}^{0}}T}}} \right. \kern-0em} {{{R}^{0}}T}})} \right\}\prod\limits_{j = 1}^N {{{{\left( {\frac{{{{X}_{j}}p}}{{{{R}^{0}}T}}} \right)}}^{{\nu _{{j,k}}^{ + }}}}} ,\,\,\,\,i = 1,...,N, \\ \end{gathered} $

где $\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],$
(16)

где $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.

Рис. 1.

Схема расчетной области.

Мгновенное положение фронта пламени в момент времени 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{:}}$

$u = {{\left[ {{{{(u_{x}^{'})}}^{2}} + {{{(u_{y}^{'})}}^{2}}} \right]}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}.$

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

(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, 1823], а также с расчетами по программе [24] для более широкого интервала объемного содержания водорода: от 6% до 75%. Отметим, что в работе [18] представлены экспериментальные данные как по ламинарным, так и по турбулентным пламенам, полученные для смесей с объемным содержанием Н2 от 9.09% до 23.09%. Ниже именно эти данные [18] используются для сравнения расчетных и измеренных скоростей турбулентного горения смесей разного состава от интенсивности турбулентности. Для водородно-воздушных смесей с содержанием водорода Н2 < 9.0% получить cтационарное решение для плоского ламинарного пламени не удалось.

Рис. 2.

Сравнение расчетных и измеренных значений скорости распространения ламинарного пламени un от состава водородно-воздушной смеси. Начальная температура Т0 = 293 К, давление атмосферное. Обозначения: 1 – расчеты данной работы, опытные точки 2 – [3], 3 – [18], 4 – [19], 5 – [20], 6 – [21], 7 – [22], 8 – [23], кривая – расчеты [24].

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

Рис. 3.

Расчетные положения фронта турбулентного пламени в разные моменты времени (t = 0–0.8 мс). Водородно-воздушная смесь, коэффициент избытка горючего Ф = 0.5 (объемное содержание Н2 – 17.36%), начальная температура Т0 = 293 К, давление атмосферное, характеристики турбулентности: $\bar {u}$ = 130 см/с, L' = 1 см, L" = 0.7 см, τ = 0.010 с.

Несмотря на относительно короткое физическое время расчета (t < τ), полученное значение скорости распространения усредненного турбулентного пламени почти постоянно во времени (рис. 4). По положениям фронта пламени, усредненным по нескольким реализациям (обычно 5–7 реализаций), можно определить видимую скорость его распространения U при данной пульсационной скорости. Если учесть тепловое расширение продуктов реакции, можно определить скорость распространения турбулентного пламени:

${{u}_{t}} = U\frac{{{{T}_{0}}}}{{{{T}_{1}}}}\frac{{{{m}_{0}}}}{{{{m}_{1}}}},$
Рис. 4.

Расчетная зависимость расстояния, пройденного усредненным фронтом турбулентного пламени, от времени. Водородно-воздушная смесь, коэффициент избытка горючего Ф = 0.5, начальная температура Т0 = 293 К, давление атмосферное, характеристики турбулентности: $\bar {u}$ = = 660 см/с, L' = 1 см, L" = 0.7 см, τ = 0.010 с.

где 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.$ Именно это условие позволяет сравнивать результаты двумерных расчетов с “трехмерным” экспериментом. Из приведенного сравнения расчетных и измеренных результатов следует вывод об их удовлетворительном согласии.

Рис. 5.

Зависимости скорости распространения турбулентного пламени ut от среднеквадратичной пульсационной скорости $\bar {u}.$ Черные кружки – расчет, светлые кружки с линиями – эксперимент [18], водородно-воздушные смеси с содержанием Н2 = 23.09% (а), 13.51% (б), 17.36% (в) и 9.09% (г), начальная температура Т0 = 293 К, давление атмосферное, характеристики турбулентности: L' = 1 см, L" = 0.7 см, τ = 0.010 с.

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

Рис. 6.

Расчетные профили температуры в ламинарном (1) и турбулентном (2) пламенах. Водородно-воздушная смесь, коэффициент избытка горючего Ф = 0.5, начальная температура Т0 = 293 К, давление атмосферное. Для турбулентного пламени характеристики турбулентности: $\bar {u}$ = 1100 см/с, L' = 1 см, L" = = 0.7 см, τ = 0.010 с.

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

Рис. 7.

Расчетные профили температуры и концентраций некоторых веществ в турбулентном пламени. Водородно-воздушная смесь, коэффициент избытка горючего Ф = 0.5, начальная температура Т0 = 293 К, давление атмосферное, характеристики турбулентности: $\bar {u}$ = 1100 см/с, L' = 1 см, L" = 0.7 см, τ = 0.010 с.

Рис. 8.

Расчетные зависимости концентрации гидроксила от среднеквадратичной пульсационной скорости $\bar {u}.$ Водородно-воздушные смеси разного состава, начальная температура Т0 = 293 К, давление атмосферное, характеристики турбулентности: L' = 1 см, L" = 0.7 см, τ = 0.010 с.

Хорошо известно, что с повышением давления скорость распространения ламинарного пламени уменьшается, а скорость распространения турбулентного пламени либо остается постоянной, либо несколько увеличивается [5]. Последнее обнаруживается и в наших расчетах распространения турбулентного пламени при разных давлениях. На рис. 9 сравниваются расчетные зависимости для скоростей распространения ламинарного и турбулентного пламен.

Рис. 9.

Расчетные зависимости скоростей распространения ламинарного, un, и турбулентного, ut, пламен от давления. Водородно-воздушная смесь, Ф = 0.5, Т0 = = 293 К, характеристики турбулентности (для ut): $\bar {u}$ = = 510 см/с, L' = 1 см, L" = 0.7 см, τ = 0.010 с.

В заключение отметим один результат, который важен для понимания структуры зоны горения в турбулентном пламени: во всех проведенных расчетах поверхность фронта пламени была односвязной. Даже при максимальных интенсивностях турбулентности расчеты не показали наличия “островков” свежей горючей смеси, окруженных продуктами горения, как это было описано в [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)$ и разыгрывается значение пульсационной скорости uxx1, 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) вычисляется

$\left\langle {{{u}_{у }}\left( {0,1,0} \right)} \right\rangle = {{u}_{у }}\left( {0,0,0} \right) + \left\{ {\left[ {{{\partial {{u}_{у }}} \mathord{\left/ {\vphantom {{\partial {{u}_{у }}} {\partial у }}} \right. \kern-0em} {\partial у }}} \right]\left( {0,0,0} \right)} \right\}\Delta {{у }_{1}};$

– в точке (х = 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)〉, определяется корректирующий коэффициент

${{a}_{u}} = {{\left\langle {{{u}_{у }}(N,{\text{1}},0)} \right\rangle } \mathord{\left/ {\vphantom {{\left\langle {{{u}_{у }}(N,{\text{1}},0)} \right\rangle } {\bar {u}}}} \right. \kern-0em} {\bar {u}}},$

находятся скорректированные значения 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) вычисляется

$\left\langle {{{u}_{х }}\left( {1,1,0} \right)} \right\rangle = {{u}_{х }}\left( {0,1,0} \right) + \left\{ {\left[ {{{\partial {{u}_{х }}} \mathord{\left/ {\vphantom {{\partial {{u}_{х }}} {\partial х }}} \right. \kern-0em} {\partial х }}} \right]\left( {1,1,0} \right)} \right\}\Delta {{х }_{1}};$

– в точке (х = 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)〉, определяется корректирующий коэффициент

${{a}_{{ux}}} = {{\left\langle {{{u}_{х }}(N,{\text{1}},{\text{ }}0)} \right\rangle } \mathord{\left/ {\vphantom {{\left\langle {{{u}_{х }}(N,{\text{1}},{\text{ }}0)} \right\rangle } {\bar {u}}}} \right. \kern-0em} {\bar {u}}},$

находятся скорректированные значения 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) пульсаций скорости:

$\rho {{u}_{k}} = {{\rho }_{0}}{{u}_{{k,0}}},$

где uk, 0 – проекции мгновенной пульсационной скорости на ось хk при начальной плотности. Этим уравнением определяется изменение интенсивности турбулентности в процессе горения, а разыгрывание пульсаций скорости требуется только при начальной температуре.

Согласно алгоритму моделирования синтетического поля турбулентности значение среднеквадратичной пульсационной скорости $\bar {u}$ постоянно во времени, а для всего расчетного пространства постоянны характеристики турбулентности L', L" и τ.

Список литературы

  1. Damkoler G. // Z. Elektrochem. 1940. B. 46. S. 601.

  2. Щелкин К.И. // ЖТФ. 1943. Т. 13. С. 520.

  3. Lewis B., Elbe G. Combustion, Flames and Explosions of Gases. Orlando: Acad. Press. 1987.

  4. Соколик А.С. Самовоспламенение, пламя и детонация в газах. М.: Изд-во АН СССР, 1960.

  5. Щетинков Е.С. Физика горения газов. М.: Наука, 1965.

  6. Bernard J.A., Bradley J.N. Flame and Combustion. L.–N.Y.: Chapman. Hall, 1985.

  7. Echekki T., Chen J.H. // Combust and Flame. 2003. V. 134. P. 169.

  8. Bell J.B., Day M.S., Grcar J.F. // Proc. Combust. Inst. 2002. V. 29. P. 1987.

  9. Bell J.B., Cheng R.K., Day M.S., Shepherd I.G. // Ibid. 2006. V. 31. P. 1309.

  10. Aspden A.J., Day M.S., Bell J.B. // Combust and Flame. 2016. V. 166. P. 266.

  11. Басевич В.Я., Володин В.П., Когарко С.М., Перегудов Н.И. // Хим. физика. 1982. Т. 1. № 8. С. 1130.

  12. Басевич В.Я., Беляев А.А., Фролов С.М., Басара Б. // Горение и взрыв. 2017. Т. 10. № 1. С. 4.

  13. Вильямс Ф.А. Теория горения. М.: Наука, 1971.

  14. Годунов С.К., Рябенький В.С. Разностные схемы. М.: Наука, 1977.

  15. Басевич В.Я., Беляев А.А., Посвянский В.С., Фролов С.М. // Хим. физика. 2013. Т. 32. № 4. С. 87.

  16. 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.

  17. Рид Р., Праусниц Дж., Шервуд Т. Свойства газов и жидкостей. Л.: Химия, 1982.

  18. Карпов В.П., Северин Е.С. // Физика горения и взрыва. 1980. Т. 16. № 1. С. 45.

  19. Козаченко Л.С. Дис. … д-ра физ.-мат. наук. М.: ИХФ АН СССР, 1954.

  20. Manton J., Milliken B.B. // Proc. Gas Dynamics Symp. (Aerothermochem.) Northwestern Univ., 1956. P. 151.

  21. Andrews G.E., Bradley D. // Combust and Flame. 1973. V. 20. P. 77.

  22. Iijima T., Takeno T. // Ibid. 1986. V. 65. P. 35.

  23. Dowdy D.R., Smith D.B., Taylor S.C., Williams A. // Proc. Combust. Inst. 1990. V. 23. P. 325.

  24. Беляев А.А., Посвянский В.С. // Алгоритмы и программы. Информ. бюлл. гос. фонда алгоритмов и программ СССР. 1985. № 3. С. 35.

  25. Басевич В.Я., Когарко С.М. // Физика горения и взрыва. 1985. Т. 21. № 5. С. 12.

Дополнительные материалы отсутствуют.