Теплофизика высоких температур, 2019, T. 57, № 1, стр. 106-113

ВВЕДЕНИЕ

При проектировании различных технических устройств широко применяются инженерные расчеты на основе решения стационарных осредненных уравнений Рейнольдса (Reynolds Averaged Navier-Stokes equations – RANS). В зависимости от особенностей течения расчеты методом RANS с большей или меньшей точностью позволяют оценить необходимые параметры потока. В [1] отмечено, что наихудшая точность расчетов наблюдается для течений с обширными отрывными зонами в случае отсутствия горения в потоке. Для повышения точности были выработаны следующие рекомендации [2]: 1) использовать нестационарную систему уравнений Рейнольдса (Unsteady RANS–URANS); 2) трактовать начальные условия на входной границе расчетной области как нестационарные и имеющие случайный характер; 3) вводить специальный дополнительный источниковый член в модель турбулентности, используемую для замыкания системы уравнений Рейнольдса; 4) применять центральные разности при аппроксимации производных с целью снижения устойчивости системы разностных уравнений; 5) проводить статистическое осреднение полученного нестационарного решения.

Перечисленные рекомендации реализованы, например, в программе Fluent из пакета прикладных программ ANSYS-CFD [2] только для расчета трехмерных течений газа, в том числе с горением. Однако расчеты течений с горением в трехмерной постановке с учетом детальных кинетических механизмов сопряжены с высокими вычислительными затратами, что сильно затрудняет проведение исследований. Для инженерных приложений более реальным представляется частный случай двумерного моделирования, но он не включен в “меню” программы Fluent.

Данная работа посвящена разработке методики моделирования турбулентного горения в двумерной постановке в рамках программы Fluent и демонстрации применения этой методики для расчета V-образного пламени.

МЕТОДИКА ДВУМЕРНОГО МОДЕЛИРОВАНИЯ ТУРБУЛЕНТНОГО ГОРЕНИЯ МЕТАНОВОЗДУШНОЙ СМЕСИ

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

Методика двумерного моделирования турбулентного горения основана на URANS-подходе. В двумерном случае этот подход не поддерживается ни одной из двух имеющихся в программе Fluent моделей случайной генерации искусственных пульсаций скорости на входной границе расчетной области (эти модели доступны только в трехмерном случае). Поэтому в данной работе используется другая модель [3], в которой нестационарные профили компонент скорости на входе в расчетную область определяются при помощи алгоритмов создания искусственной турбулентности. Данные алгоритмы позволяют задавать в набегающем потоке необходимые значения как интенсивности пульсаций, так и интегрального масштаба турбулентности. В техническом плане для этого применяется предусмотренная в программе Fluent возможность использовать UDF-функции (User Defined Function).

Заметим, что задание скорости на входной границе расчетной области (так называемые “жесткие” граничные условия, по терминологии Fluent “velocity-inlet”) при URANS-подходе к дозвуковому течению допустимо (в рамках программы Fluent) при условии использования особой зависимости для плотности. Вместо уравнения состояния для идеального газа (“ideal-gas”) в этом случае используется другая зависимость “incompressible-ideal-gas” (подробнее см. [2]).

В качестве модели турбулентности используется широко распространенная разновидность k–ε-модели – “Realizable” [2]. Никаких специальных дополнительных источниковых членов эта модель не содержит, в отличие от упомянутой во Введении рекомендации 3.

Остановимся подробнее на задании граничных условий. Ранее опубликованные расчетные результаты [4] показали, что учет турбулентности набегающего потока только посредством граничных условий (стационарных) в модели турбулентности позволяет правильно предсказать угол наклона средней линии фронта пламени, но не позволяет правильно предсказать расширение фронта по длине. При этом учет турбулентности набегающего потока только посредством нестационарных граничных условий для скорости позволяет правильно предсказать расширение фронта по длине [3]. Но в этом случае приходится решать уравнения Навье–Стокса на чрезвычайно мелкой сетке, что, как известно, сопряжено с большими вычислительными трудностями. В данной работе рассматривается вариант, когда учет турбулентности набегающего потока осуществляется посредством и граничных условий (стационарных) в модели турбулентности, и нестационарных граничных условий для скорости. Как будет показано ниже, результаты расчетов в этом случае с хорошей точностью описывают и угол наклона фронта пламени, и его утолщение по мере удаления от места стабилизации пламени. Возможным способом устранения некорректности, связанной с двойным учетом турбулентности, может быть вариант, при котором часть энергии турбулентности набегающего потока учитывалась бы в граничных условиях для модели турбулентности, а оставшаяся часть – в граничных условиях для скорости (в пропорции, например, 50/50).

В качестве модели горения используется “квазиламинарная” (QL) модель [5], в которой применяется модель турбулентности для описания эффектов смешения, но не учитывается взаимодействие химических источников с турбулентностью. В программе Fluent этой модели горения соответствует опция “Laminar Finite-Rate” в окне “Turbulence-Chemistry Interaction”. Заметим, что в случае горения метановоздушной смеси расчеты модельных камер сгорания методом RANS с использованием QL-модели горения удовлетворительно описали эмиссию как оксидов азота, так и оксида углерода [4, 6]. Но, кроме эмиссионных характеристик, важными являются и характеристики фронта горения: положение середины фронта и его толщина. Расчеты V-образного пламени методом RANS с использованием QL-модели горения удовлетворительно описывают положение середины фронта горения, где уровень температуры соответствует 50% от ее изменения на всем фронте [4]. Что же касается толщины фронта пламени, то в эксперименте она линейно увеличивается по мере удаления от места воспламенения [7], а расчеты этот эффект, к сожалению, не описывают [4]. В этом состоит основной недостаток использования метода RANS с QL-моделью горения.

В качестве кинетического механизма используется редуцированный механизм DRM-19 [8], включающий 84 реакции для 21 компонента, причем два из них – Ar и N₂ – являются инертными разбавителями, не участвующими в химических реакциях. Данный механизм представляет собой усеченную версию известного детального кинетического механизма GRI-Mech 1.2 [9], содержащего 177 реакций с участием 32 компонентов.

Из-за упомянутых выше допущений полученные в данной работе результаты не могут считаться строгими (хотя и с неплохой точностью описывают характеристики плоского V-образного пламени, полученные в эксперименте) и носят, скорее, качественный характер.

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

Рассмотрим натекание турбулентного потока гомогенной метановоздушной смеси на нагретое плохо обтекаемое тело. V-образное пламя стабилизируется на этом теле и далее распространяется от него в набегающий поток несгоревшей смеси. Схема течения, соответствующая рассматриваемой задаче, приведена на рис. 1 (U – вектор осредненной скорости).

Расчетная область представляет собой равнобедренную трапецию. Размеры ее оснований составляют AB = 15 см и DC = 11 см, высота – 11 см. Размер торца плохо обтекаемого тела, расположенного в сечении x = 0, равен 2 мм. На входе в расчетную область x = 0 (линия AB на рис. 1) при –1 < y < 1 мм задаются условия стенки, нагретой до 2500 К, а на остальном участке (y > 1 мм и y < ‒1 мм) – нестационарные профили продольной u и вертикальной v компонент скорости при помощи специальной процедуры случайных чисел, моделирующей искусственную турбулентность (подробнее см. [3]). На линии BCDA, ограничивающей расчетную область, задаются условия для выходной границы.

Параметры набегающей на фронт горения смеси соответствуют условиям эксперимента [7]. Скорость бедной метановоздушной смеси составляет U_∞ = 3 м/с, давление P = 1 атм, температура T_∞ = 298 К, коэффициент избытка воздуха α = 1.43, интегральный масштаб турбулентности L_t = 3.5 мм, интенсивность турбулентных пульсаций скорости в продольном направлении I = 7%, в поперечном I = 5.5%. В расчете задавалось среднее значение I = 6.25%.

Численное моделирование проводилось в программе Fluent [10] из пакета прикладных программ ANSYS-CFD. Вязкость задавалась с учетом многокомпонентности среды (ideal-gas-mixing-law). Расчеты выполнялись со вторым порядком точности по времени и третьим по пространству. Использовалась одноблочная структурированная сетка. Число ячеек составляло N = 16 500, а размер ячейки ∼1 мм. Заметим, что использовать центральные разности в соответствии с рекомендацией 4 не представляется возможным, так как в рассматриваемом двумерном случае эта опция отсутствует.

В процессе нестационарного расчета осуществлялся сбор статистики в течение некоторого времени τ. Это время соотносилось с некоторым характерным временем t *, в качестве которого взято время пролета расчетной области вмороженной в набегающий поток частицей при отсутствии горения. В расчетах сбор статистики осуществлялся в течение примерно τ = 15t *.

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

Сначала рассмотрим режим, соответствующий экспериментальным условиям [7]. На рис. 2 приведено рассчитанное поле статистически осредненной температуры. Проведем обработку поля температур с помощью следующей процедуры. Перепад температуры во фронте составил 1522 К: от температуры набегающего потока T_∞ = = 298 К до максимальной температуры за фронтом горения T_max = 1820 К. Рассмотрим три изолинии температуры, характеризующие холодную границу (10% от перепада температуры), среднюю линию (50%) и горячую границу (90%). Расположение этих трех изолиний на расчетной плоскости приведено на рис. 3а. Для сравнения на рис. 3б показано положение этих изолиний при аналогичном расчете, но без задания искусственной турбулентности на входной границе. Видно, что в случае задания искусственной турбулентности на входе ширина фронта горения линейно увеличивается по мере удаления от места стабилизации пламени, а в случае отсутствия пульсаций скорости ширина фронта горения примерно постоянна.

На рис. 4 представлено сравнение рассчитанных значений толщины фронта горения δ, которая определяется как расстояние между холодной и горячей границами фронта пламени, с данными эксперимента [7]. Видно, что наблюдается удовлетворительное соответствие рассчитанного и измеренного темпов роста толщины фронта горения. Отличие в абсолютном значении δ при малых значениях x обусловлено наличием в расчете плохо обтекаемого тела с донным торцом толщиной 2 мм.

Оценим темп утолщения фронта горения. Для этого определим углы наклона расчетных границ β к направлению набегающего потока. Для оценок используем положения границ при x = 6 и 11 см. Получены следующие значения:

холодная граница – β_cold = 13.5° для y > 0 и β_cold = = 13.3° для y < 0,

средняя линия – β_mid = 11.07° для y > 0 и β_mid = = 11.55° для y < 0,

горячая граница – β_hot = 9.61° для y > 0 и β_hot = = 8.65° для y < 0.

Тогда угол утолщения фронта горения Δβ = = β_cold – β_hot = 3.9° для y > 0 и Δβ = 4.65° для y < 0, т.е. расчетные значения Δβ для верхней и нижней полуплоскостей несколько различаются. Будем использовать среднеарифметические значения углов утолщения фронта 〈Δβ〉 = 4.27° и положения средней линии фронта 〈β_mid〉 = 11.3°. В эксперименте [7] 〈β_mid〉 = 12.2°. Это значение отличается от расчетного на 8%. Экспериментальное значение угла утолщения фронта определялось путем обработки экспериментальных данных для зависимости толщины фронта от расстояния, представленных в [7] и на рис. 4 данной работы. Это значение, рассчитанное по данным о ширине фронта пламени при x = 2 и 10 см, составило 〈Δβ〉 = 4.4°, что отличается от расчетного значения на 3.4%. Таким образом, можно констатировать удовлетворительное соответствие расчетных значений экспериментальным данным как по положению фронта горения, так и по его утолщению.

Турбулентная скорость горения S_turb в случае V-образного пламени определяется следующим образом. В точке пересечения вертикали x = 8.5 см с холодной границей фронта горения находим вектор скорости потока и проецируем его на перпендикуляр к средней линии фронта горения. Эту проекцию и будем полагать турбулентной скоростью горения. Полученные таким образом значения S_turb для рассматриваемого случая равны 0.192 м/с для y > 0 и 0.233 м/с для y < 0. Среднеарифметическое значение турбулентной скорости горения равно 0.21 м/с. Ламинарная скорость горения в рассматриваемом случае составляет S_lam = 0.19 м/с [7]. Таким образом, отношение турбулентной скорости горения к ламинарной равно 〈S_turb〉/S_lam = 1.12.

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

Влияние интенсивности турбулентных пульсаций скорости набегающего потока. Кроме рассмотренного выше режима с I = 6.25%, были проведены расчеты для I = 10 и 14%. Границы фронтов горения для различных значений I представлены на рис. 5. Видно, что увеличение интенсивности турбулентных пульсаций скорости I приводит к более интенсивному утолщению фронта горения. Кроме того, при I = 14% (рис. 5в) наблюдается существенное смещение положения горячей границы фронта пламени вниз по потоку по сравнению со случаями с I = 6.25 и 10% (рис. 5а и 5б). Анализ расчета течения с I = 14% выявил факты временного срыва горения с его последующим восстановлением. Один из таких моментов представлен на рис. 6. Видимо, этим обстоятельством и объясняется зафиксированная аномалия в форме горячей границы фронта при I = 14%.

Рассчитанные значения скорости турбулентного горения, положения и толщины фронта при различных I представлены в табл. 1. Увеличение интенсивности пульсаций скорости приводит к росту турбулентной скорости горения. На рис. 7 показана рассчитанная зависимость S_turb/S_lam от среднеквадратичной пульсации скорости u'/S_lam. В настоящее время существует несколько эмпирических зависимостей, учитывающих влияние интенсивности пульсаций на турбулентную скорость горения (см., например, [11]). Одна из них – модифицированная в [11] эмпирическая зависимость [12] – имеет следующий вид:

где Re_t = u'L_t/ν – турбулентное число Рейнольдса, u' – среднеквадратичная пульсация скорости, ν – кинематическая вязкость. В расчете величина ν для метановоздушной смеси с коэффициентом избытка топлива ϕ = 0.7 и T_∞ = 298 К принималась равной 1.6 × 10^–5 м²/с. Видно (рис. 7), что расчетные данные качественно согласуются с эмпирической зависимостью (1). Этот факт можно расценивать как еще одно косвенное подтверждение применимости рассматриваемого подхода к моделированию турбулентного горения топливовоздушной смеси.

Увеличение интенсивности турбулентных пульсаций скорости приводит также к увеличению угла средней линии фронта. При увеличении I от 6.25 до 14% (т.е. на 124%) 〈β_mid〉 возрастает с 11.3° до 14.6° (30%). Увеличение I приводит к немонотонному изменению ширины фронта. Если при росте I от 6.25 до 10% значение 〈Δβ〉 увеличивается от 4.27° до 6.76° (58%), то дальнейшее увеличение интенсивности от 10% до 14% ведет уже к уменьшению 〈Δβ〉 от 6.76° до 4.7° (30%). Скорее всего, этот немонотонный характер зависимости обусловлен отмеченным выше эффектом срыва горения с последующим его возобновлением.

Влияние интегрального масштаба турбулентности набегающего потока. Дополнительно к упомянутому выше режиму с L_t = 3.5 мм рассмотрим режимы с L_t = 2 и 5 мм. Интенсивность турбулентных пульсаций в этих расчетах полагалась равной I = 6.25%. Результаты расчетов представлены в табл. 2. Видно, что увеличение L_t набегающего потока в 2.5 раза приводит к совсем незначительному уменьшению относительной скорости турбулентного горения от 1.126 до 1.107, т.е. всего на 1.7%. Если воспользоваться эмпирической зависимостью (1), то изменение L_t от 2 до 5 мм ведет к уменьшению 〈S_turb〉/S_lam на 2%. Расчетное значение (1.7%) практически совпадает с эмпирической оценкой (2%). Заметим, что увеличение I в сопоставимые 2.24 раза (от 6.25 до 14%) привело к гораздо более сильному изменению 〈S_turb〉/S_lam, а именно к увеличению на 36% (см. табл. 1). Столь же малое влияние L_t отмечается и на положении фронта горения. Так, угол средней линии фронта 〈β_mid〉 изменился лишь на 1.3% (см. табл. 2). Однако изменение L_t сильно повлияло на ширину фронта. При увеличении L_t в 2.5 раза, 〈Δβ〉 выросло в 1.8 раза.

Влияние неоднородности состава топливовоздушной смеси. Рассмотрим ситуацию, когда набегающий поток топливовоздушной смеси не полностью перемешан, т.е., кроме флуктуаций скорости, в нем присутствуют неоднородности состава. Предположим по аналогии со скоростью потока, что неоднородность состава топливовоздушной смеси характеризуется интенсивностью ${{I}_{{{\text{C}}{{{\text{H}}}_{4}}}}}$ и интегральным масштабом пульсаций мольной доли метана ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}}.$ Кроме уже рассмотренного случая однородной топливовоздушной смеси (${{I}_{{{\text{C}}{{{\text{H}}}_{4}}}}}$ = 0), рассчитаны варианты с не полностью перемешанной смесью с ${{I}_{{{\text{C}}{{{\text{H}}}_{4}}}}}$ = 2, 4 и 6%. Интегральный масштаб пульсаций мольной доли метана ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}}$ составлял в этих расчетах 3.5 мм. Интегральный масштаб и интенсивность турбулентных пульсаций скорости задавались такими же, как в эксперименте [7]: L_t = 3.5 мм и I = 6.25%. Результаты расчетов представлены в табл. 3. Видно, что относительная скорость турбулентного горения неоднородной смеси 〈S_turb〉/S_lam отличается от относительной скорости горения однородной смеси не более чем на 2.6%. Угол средней линии фронта 〈β_mid〉 в случае неоднородной смеси увеличивается не более чем на 1% по сравнению с однородной смесью, а угол уширения фронта горения 〈Δβ〉 уменьшается не более чем на 20%.

Для случая с ${{I}_{{{\text{C}}{{{\text{H}}}_{4}}}}}$ = 6% рассмотрены режимы с различными значениями ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}}$ (2, 3.5 и 5 мм). Результаты расчетов представлены в табл. 4. Видно, что с увеличением ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}}$ скорость турбулентного горения неоднородной смеси 〈S_turb〉/S_lam уменьшается. В целом эта характеристика изменилась не более чем на 10%. C увеличением ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}}$ также уменьшается 〈β_mid〉. Изменение составило не более 7%. Угол уширения фронта горения 〈Δβ〉 увеличивается с ростом ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}},$ но не более чем на 3.3%. При этом интегральный масштаб пульсаций мольной доли метана ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}}$ в этих расчетах изменялся в 2.5 раза.

Таким образом, можно констатировать, что все рассмотренные характеристики горения мало изменяются при переходе от однородной к не полностью перемешанной смеси с интенсивностью пульсаций мольной доли топлива ${{I}_{{{\text{C}}{{{\text{H}}}_{4}}}}}$ по крайней мере до 6% и интегральным масштабом пульсаций мольной доли топлива ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}}$ до 5 мм. С практической точки зрения это означает, что горение полностью и не полностью перемешанной смесей происходит практически одинаково.

Ранее проведенные теоретические, численные и экспериментальные исследования имеют противоречивый характер. Например, экспериментальные исследования [13] свидетельствуют о более чем трехкратном увеличении скорости турбулентного горения пропановоздушной смеси при увеличении степени неоднородности смеси (см. рис. 10 в [13]). В этом эксперименте осредненный коэффициент избытка топлива равен ϕ = 0.7. Полученные в настоящей работе результаты не подтверждают этот эффект. В экспериментальной работе [14] не зафиксировано увеличения турбулентной скорости горения вследствие неоднородности метановоздушной смеси. DNS-расчеты (Direct Numerical Simulation), выполненные в [15], также не выявили увеличения турбулентной скорости горения вследствие случайных флуктуаций состава стехиометрической смеси (ϕ = 1), хотя в случае бедной смеси (ϕ = 0.7) такое увеличение имело место.

Влияние конфигурации расчетной области (эффект горизонтально ориентированной разделительной пластины). Пластина размещается за плохо обтекаемым телом вдоль оси x. Пограничный слой, образующийся на пластине, не учитывается. В качестве граничного условия используется условие симметрии. Остальные параметры течения соответствуют тем, что использовались при моделировании эксперимента [7]. Расположение на расчетной плоскости холодной, горячей и средней температурных границ в этом случае приведено на рис. 8, а рассчитанные значения 〈S_turb〉/S_lam, положения и толщины фронта представлены в табл. 5. Из сравнения видно, что наличие пластины способствует подавлению пульсаций скорости, а это приводит к уменьшению 〈S_turb〉/S_lam на 12% и к снижению 〈β_mid〉 на 8%. При этом 〈Δβ〉 уменьшается значительно – почти в три раза.

ЗАКЛЮЧЕНИЕ

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

Проведенный на основе двумерных нестационарных расчетов параметрический анализ влияния различных параметров набегающего потока на характеристики V-образного пламени показал:

1) качественное соответствие расчетных данных эмпирической зависимости [12] по влиянию интенсивности турбулентных пульсаций скорости на турбулентную скорость горения;

2) увеличение турбулентной скорости горения и угла средней линии фронта пламени с ростом интенсивности пульсаций скорости (при увеличении интенсивности пульсаций скорости от 6.25 до 14% скорость турбулентного горения увеличивается на 36%, а угол средней линии фронта пламени на 30%);

3) незначительное влияние интегрального масштаба турбулентности на скорость турбулентного горения (при увеличении интегрального масштаба турбулентности в 2.5 раза скорость турбулентного горения уменьшается не более чем на 2%);

4) значительное влияние интегрального масштаба турбулентности на изменение угла уширения фронта (при увеличении интегрального масштаба турбулентности в 2.5 раза угол уширения фронта увеличивается в 1.8 раза);

5) незначительное влияние полноты смешения метана с воздухом (с интенсивностью пульсаций мольной доли топлива ${{I}_{{{\text{C}}{{{\text{H}}}_{4}}}}}$ по крайней мере до 6% и интегральным масштабом пульсаций мольной доли топлива ${{L}_{{t,{\text{C}}{{{\text{H}}}_{4}}}}}$ до 5 мм) на характеристики турбулентного горения;

6) уменьшение как турбулентной скорости горения (на 12%), так и угла уширения фронта (в три раза) при введении в поток разделительной пластины, ориентированной по набегающему потоку.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 17-08-01423).

Авторы благодарят П.Д. Токталиева за содействие в обеспечении вычислительными ресурсами.