CC BY
94
УДК 004.942+621.452.3
МОДЕЛИРОВАНИЕ СТРУКТУРЫ ПОТОКА ЗА СТАБИЛИЗАТОРОМ ПЛАМЕНИ МЕТОДОМ КРУПНЫХ ВИХРЕЙ
© 2013 С.Г. Матвеев, И.А. Зубрилин
Самарский государственный аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет)
Поступила в редакцию 02.12.2013
В статье представлены и проанализированы результаты расчётов потока за стабилизатором пламени при использовании подходов осреднения уравнений Навье-Стокса по Рейнольдсу и метода моделирования крупных вихрей. Показано, что подход моделирования турбулентности, основанный на осреднении по Рейнольдсу не способен прогнозировать крупные нестационарные вихревые структуры. Для подобного типа течений перспективным является использование метода моделирования крупных вихрей, который показал более близкие к экспериментальным данным результаты. Исследована возможность использования метода моделирования крупных вихрей в сочетании с глобальными и редуцированными механизмами химических реакций для расчёта бедного срыва пламени пропано-воздушной смеси.
Ключевые слова: стабилизатор пламени, вычислительная газовая динамика, модели турбулентности, моделирование крупных вихрей.
Устойчивое горение пламени во всем диапазоне рабочих режимов двигателя является одним из основных требований к камерам сгорания (КС) ГТД [1, 2]. Обеспечение устойчивого горения становится особенно актуальным при использовании технологии сжигания заранее подготовленных бедных смесей на переходных режимах работы двигателя, при которых коэффициент избытка воздуха а на выходе из КС может достигать значений я=15...20 [2]:
<=tV
где СЕ - расход воздуха, Ст - расход топлива,£0-стехиометрический коэффициент.
Экспериментальное исследование пределов устойчивой работы на переходных режимах работы двигателя связано с высокими материальными и временными затратами. Использование методов вычислительной газовой динамики на этапе проектирования позволяет сузить диапазон наиболее оптимальных конструктивных решений и сократить объем доводочных испытаний.
Наиболее распространенными способами моделирования турбулентных течений, хорошо зарекомендовавшими себя в последние несколько десятилетий является решение осредненных по Рейнольдсу уравнений Навье-Стокса (Reynolds Average Navier-Stokes, RANS) и метода
Матвеев Сергей Геннадьевич, кандидат технических наук,
доцент кафедры теплотехники и тепловых двигателей.
E-mail: [email protected]
Зубрилин Иван Александрович,
теплотехники и тепловых двигателей.
Email: [email protected]
моделирования крупных вихрей (Large Eddy Simulation, LES) [3, 4, 5]. Метод RANS основан на тотальном осреднении всего спектра турбулентных пульсаций и его представлении в виде характеристик турбулентности, таких как, кинетическая энергия турбулентности и скорость её диссипации. Его преимуществами являются относительно невысокие требования к пространственной дискретизации и параметрам вычислительной техники. Недостатками - отсутствие возможности моделирования вихревых структур с размерами порядка макромасштаба потока, наличие эмпирических коэффициентов при замыкании уравнений и, как следствие, отсутствие приемлемой универсальности. Метод крупных вихрей основан на математическом разделении потока на крупновихревое и мелковихревое движение с помощью процедуры фильтрации. Влияние мелковихревого движения описывается приближенными, так называемыми подсеточными моделями. В качестве фильтра в основном выступает локальный размер конечно-элементной сетки. Преимуществами метода крупных вихрей являются возможность получения параметров среднего течения, а также спектральных характеристик, пространственных и временных масштабов турбулентности. Недостатком подхода LES являются более высокие требования к вычислительной технике.
Целью данной работы является сравнение методов моделирования турбулентности на основе решения осредненных по Рейнольдсу уравнений Навье-Стокса в нестационарной постановке (Unsteady Reynolds Average Navier-Stokes, URANS)
и метода моделирования крупных вихрей (Large Eddy Simulation, LES) при исследовании течения за уголковым стабилизатором пламени.
В качестве предмета исследования выбрано плохообтекаемое тело в виде правильной треугольной призмы с длиной ребра основания h = 25 мм, установленное в канале квадратного сечения 50х50 мм. Скорость набегающего потока равна 10 м/с, профиль скорости - равномерный. Интенсивность турбулентности на входе - 2%. Схема экспериментального участка представлена на рис. 1. Исходные данные и результаты экспериментальных измерений взяты из статьи [6]. В расчёте длина канала до стабилизатора принята равной 4h и вниз по потоку от стабилизатора 10h.
Расчёт проводился в трехмерной нестационарной постановке в программном комплексе ANSYS Fluent 14.5. Для создания блочно-струк-турированной конечно-элементной сетки геометрическая модель расчётной области была разделена на 7 блоков (рис. 2).
На первом этапе проводился расчёт изотермического течения воздух. При расчёте в постановке URANS использовались следующие модели турбулентности, реализованные в ANSYS Fluent: k-e Standard, k-e Realizable, k-ю SST,
Reynolds Stress Model (RSM). При использовании URANS необходимо исследовать зависимость решения от плотности дискретизации расчётной области и размера временного шага [7]. Для этого проводятся расчёты на нескольких вариантах сеток с последовательным увеличением плотности в 2 и более раз. В нашем случае расчёт проводился для трех различных вариантов сеток, количество конечных элементов в которых равнялось 0,70 млн., 1,43 млн. и 3,69 млн. штук. Для определения зависимости решения от шага по времени использовались различные значения временных шагов: 5 мс, 1 мс, 0,1 мс.
Расчёт в постановке LES проводился на сетке 3,69 млн. с шагом по времени 0,1 мс. Исходя из формулировки подхода LES, необходимо, чтобы размер конечно элементной сетки позволял разрешить самые крупные вихревые структуры, находящиеся вне инерционного интервала спектра кинетической энергии турбулентности (рис. 3) [4]. Для развитого турбулентного течения этот участок описывается законом Колмогорова-Обухова (закон «5/3») [5]:
Е(к) = С^к"5'3
где Е - энергетический спектр, С^=1,5 - постоян-
Рис. 1. Схема экспериментального участка [6]
Рис. 2. Геометрическая модель расчётной области
1о8(к)
Рис. 3. Спектр турбулентной кинетической энергии
ная Колмогорова, к- кинетическая энергия турбулентности, Е - скорость диссипации кинетической энергии турбулентности.
На рис. 4 представлены результаты преобразования Фурье абсолютной скорости поток, по-
лученной в результате расчета изотермического течения с использованием динамической подсе-точной модели Smagorinsky-Lilly. Изменение скорости во времени записывалось для трех точек со следующими координатами: №1 3^=0,5h, J£|= 1 h; №2yI=0,5h,xa=2h; №3 уэ =0,5h, JC3=3h. Из графиков видно, что для всех точек имеется участок, соответствующий закону «5/3» Колмогорова-Обухова. Также обнаружен ярко выраженный пик на частоте 170-190 Гц, что соответствует экспериментальным данным [6].
Для выбора подхода подсеточного моделирования при использовании LES были проведены расчеты с использованием следующих моделей: динамическая модель Smagorinsky-Lilly (SL), Wall-Adapted-Local-Eddy (WALE) и WallModeled (WM) [3, 4]. На рисунке 5 представлены графики изменения осевой скорости Vx вдоль центральной линии модели (ось x) и расхода воздуха в зоне обратных токов вдоль ее длины, где G - расход через зону обратных токов (ЗОТ),
о
L О
i J.
f
Я)
L_
о.
о х
m
1,0Е+01 1,0Е+00 1,0Е-01 1,0Е-02 1,0Е 03 1,0Е-04 1,0Е-05 1,0Е-06 1,0Е-07 1,0Е-08
Л
J Пл У til ■|l ■, \-5/3
if« ш к х
И Ш
I '1 Ж
Щ
i
i
i 1
1,0Е+01 1,0Е+02 1,0Е+03
1,0Е+04 Гц
1,0Е+01
1,0Е+00
1,0Е-01
1,0Е-02
1,0Е-03
и —. s 1,0Е-04
1,0Е-05
l,0E-06
1л0Е-07
l,0E-08
Л-ЛА1А1 lAu ji \-5/3
Щ\ к 1 Х
! ' и к
i * Hl
i i Ж
i i 1
i i i II
i i
1,0Е+01 1,0Е+02 1,0Е+03
1,0Е+04 0, Гц
1,0Е+01 1,0Е+00 1,0Е-01 1,0Е-02
=г
1,0Е 03
U
^ 1,0Е-04
1,0Е-05
1,0Е-06
1,0Е-07
1,0Е-08
ljOE+Ol
л л Aito J -5/3
VsH V
V
1 I
'к.
1,0Е+02 1,0Е+03
в)
1,0Е+04 Гц
Рис. 4. Результаты быстрого преобразования Фурье для скорости в трех точках а) - точка №1, б) - точка №2, в) - точка №3
15
5
10
-5
-10
• •
'С'
f • / .V / / */ / s
ЛИ • # ^Зг ^ * Ж / у» ? А' А* ф
j
0
4
x/h
LES SL •••• LES WM
Рис. 5. Осевая скорость потока вдоль центральной линии (а) и расход воздуха в зоне обратных токов (б)
GB - расход воздуха на входе в расчётную область, 1ЗОТ - протяженность ЗОТ.
Из графика видно, что все используемые модели дают похожие результаты, и адекватно отражают осредненные по времени характеристики потока. Однако, наиболее близкие результаты к экспериментальными данным показала динамическая модель Smagorinsky-Lilly, которая и была выбрана для дальнейших расчётов и сравнения с URANS.
Одним из определяющих параметров для устройств стабилизации пламени является размер зоны обратных токов [1, 8]. Полученные размеры ЗОТ в осевом направлении для всех вариантов конечно-элементных сеток, временных шагов и моделей турбулентности представлены в табл. 1. Значения указаны в калибрах относительно длины ребра стабилизатора. Экспериментальное
значение равнялось 2,1h [6]. На рис. 6 представлены поля мгновенной скорости потока при расчёте с использованием подхода URANS и LES.
Из результатов видно, что при расчёте изотермических течений подход URANS не позволяет получить крупновихревого нестационарного течения, наблюдаемого в эксперименте. Протяженность зоны обратных токов при этом в 1,7 и более раз больше, чем при натурных исследованиях. Полученные данные согласуются с результатами, представленными в работах [3, 13, 14]. При расчёте с использованием LES размер ЗОТ отличается от эксперимента не более чем на 5%, даже на самой крупной сетке.
На следующем этапе проводилось исследование горения заранее подготовленной пропано-воздушной смеси. Расчёт процесса горения газообразного топлива описывался с использовани-
Таблица 1. Протяженность ЗОТ при расчёте на различных моделях турбулентности
Модель турбулентности Протяженность ЗОТ, h
5мс 1мс 0,1мс
Сетка £ k-e STN 3,62 3,65 3,55
№1 k-e Real 4,10 4,18 4,16
0,70 млн. <С k-w SST 4,51 4,52 4,51
В RSM 4,89 4,87 4,89
LES SL - - 2,08
Сетка оо 5 k-e STN 3,63 3,59 3,63
№2 k-e Real 4,14 4,12 4,12
1,43 млн. <с k-w SST 4,52 4,52 4,52
В RSM 5,42 4,90 4,01
LES SL - - 2,15
Сетка оо £ k-e STN 3,78 3,63 3,67
№3 k-e Real 4,08 4,12 4,10
3,69 млн. <С k-w SST 4,55 4,53 4,54
В RSM 5,20 5,19 4,90
LES SL - - 2,13
г 516+01 2.3Эе+01 2.2®е+01 2 1Зе+01 2,016+01 1 886+01 1.766+01 1 вЗе+01 1 51е+01 1 ЗЗе+01 1 .26е+01
1 13е+01 1.008+01 В.7Эе+00 7 53е+00 6.28е+00 5.026+00 3.776+00
2 516+00 1.266+00 О.ООе+ОО
L,
а)
■
322е+01 3.06е+01 2.90е+01 2.74е+01 2.58е+01 2.42е+01 2.26е+01 2.10е+01 1 93е+01 1.77е+01 1.61е+01 1 45е+01 1.29е+01 1.t3e+01 9.67е+00 8.06е+00 6.45е+00 4 84е+00 3.22е+00 1 61е+00 О.ООе+ОО
ь
б)
Рис. 6. Поля мгновенных скоростей потока: а - URANS RSM, б - LES SL
ем модели конечной скорости химическом реакции, которая предполагает расчёт скоростей реакций по уравнению Аррениуса. В данной работе использовались два механизма: глобальная реакция окисления пропана и редуцированный кинетический механизм, предложенный в работе [9], который содержит 17 компонентов и 28 реакций. Для определения границы бедного срыва пламени использовался следующий алгоритм. В начале, проводился расчёт течения изотермической заранее подготовленной смеси с б=1. Затем, за стабилизатором пламени локально увеличивалась температура потока выше необходимой для воспламенения (1500 К). Далее проводился расчёт до достижения устойчивого, периодически повторяющегося решения. Условием этого было повторяемость результатов за два последовательных промежутка времени, равных 3"ЬВП, где "Ьвп " характерное время пребывания. Согласно рекомендации [10] характерное время определяется как:
^ЕП
V
,
где V - объем расчётной области, Q расход на входе.
объемный
Результаты сравнивались по полю осреднён-ной по времени температуры и полю её среднеквадратичного отклонения. Расчёт продолжался
, затем на входе в расчётную область ступенчато увеличивался коэффициент избытка воздуха с шагом Да = ОД. Для определения момента срыва пламени на каждом временном шаге записывалась среднемассовая температура потока в трех сечениях: на расстоянии x=2h, x=5h и на выходе из расчётной области.
Изменение температуры во времени в сечении 2Ь при расчёте с использованием глобальной реакции представлено на рис. 7. Из графика видно, что бедный срыв пламени происходит при ££=2,5-2,6, экспериментальное значение а= 1,9 [6]. Это свидетельствует о том, что модель горения, основанная на глобальной реакции, завышает реакционную способность бедных смесей.
Изменение температуры во времени в трёх сечениях при расчёте с использованием кинетического механизма, предложенного в работе [9] представлено на рисунке 8. В этом случае погасание пламени в сечении х=2Ь происходит при '■ = 1,6. Падение температуры в сечениях, которые находятся вниз по потоку, происходит с неболь-
Рис. 7. Температура в сечении х=2Ь при расчёте с использованием 1 реакции
Рис. 8. Изменение температуры при расчёте с использованием редуцированной кинетической схемы [9]
шим запаздыванием. Однако, данное значение ниже, чем экспериментальное. Возможно, это связано с недостаточным пространственным разрешением, в результате чего поверхность фронта ламинарного пламени оказалось меньше, чем в реальности, как следствие занижена турбулентная скорость распространения пламени. Решением может быть использование моделей горения утолщенного фронта пламени [11], для которых при расчёте толщины фронта пламени и скорости его распространения учитывается сеточное разрешение или использование специальных процедур подсеточного моделирования [12].
Также обнаружено, что на предсрывном режиме стабилизация пламени определяется поведением крупных вихревых структур: их временем пребывания в ближнем следе за стабилизатором, частотой схода с кромок и взаимодействием между собой.
В результате можно сказать, что при использовании подхода URANS для моделирования течений с крупными трехмерными вихревыми структурами не удалось получить нестационар-
ной картины течения, протяженность ЗОТ завышена по сравнению с экспериментом в 1,7 и более раз. Более перспективным для данного типа течений является подход моделирования крупных вихрей. Данный подход показал удовлетворительное согласование с экспериментом по протяженности ЗОТ, распределению скорости вдоль центральной линии и расходу воздуха через ЗОТ. Также показана принципиальная возможность исследования явления бедного срыва пламени при заданном законе изменения б как функции от времени. Выявлено, что существенное влияние на результаты расчёта оказывает используемый кинетический механизм химических реакций. Для более точного описания поведения пламени необходимо использовать расширенные механизмы химических реакций и апробированные модели турбулентно-химического взаимодействия.
СПИСОК ЛИТЕРАТУРЫ
1. Lefebvre A.H., Ballal D. R. Gas Turbine Combustion: Alternative Fuels and Emissions. Third Edition. CRC Press, 2010. 560 p.
2. Основы конструирования авиационных двигателей и энергетических установок/ЛЛ Иноземцев, М.А. Нихамкин, В.Л. Сандрацкий. М.: Машиностроение, 2007. 396 с.
3. Волков К.Н., Емельянов В.Н. Моделирование крупных вихрей в расчетах турбулентных течений. М.: Физматлит, 2008. 368 с.
4. Снегирёв А.Ю. Высокопроизводительные вычисления в технической физике. Численное моделирование турбулентных течений : учеб. пособие. СПб.: Изд-во Политехн. ун-та, 2009. 143 с.
5. Моделирование турбулентности в расчетах сложных течений: учебное пособие / А.В. Гарбарук, М.Х. Стрелец, М.Л. Шур. СПб.: Изд-во Политехн. ун- та, 2012. 88 с.
6. Fujii, S., Fujii S. A comparison of cold and reacting flows around a bluff-body flame stabilizer // Journal of Fluids Engineering. 1981. Vol. 103, P. 803-818.
7. Verification and validation in computational fluid dynamics / L. William, W.L. Oberkampfa, T. G. Trucanob // Progress in Aerospace Sciences. 2002. Vol. 38, Issue 3, P. 209-272.
8. Талантов А.В. Горение в потоке. М.: Машиностроение, 1978. 160 с.
9. Simplified Jet-A/Air Combustion Mechanisms for Calculation of NO Emissions / K.P. Kundu, P.J. Penkot, S.L. Yang // AIAAX98-3986.
10. Menter F.R. Best Practice: Scale-Resolving Simulation in ANSYS CFD. ANSYS Inc., 2012. 70p.
11. O'Rourke P.J., Bracco F. V. Two Scaling Transformations for the Numerical Computation of Multidimensional Unsteady Laminar Flames // Journal of Computational Physics. 1979. № 33 (2). P.185 -203.
12. Poinsot T., Veynante D. Theoretical and Numerical Combustion. Third Edition. Bourdeux, France: Aquaprint, 2012. 603p.
13. Towards Modeling Lean Blow Out in Gas Turbine Flameholder Applications / W.-W. Kim, P. R. Van Slooten, R. E. Malecki, S. Syed, M. B. Colket, J. J. Lienau // J. Eng. Gas Turbines Power № 128(1) P. 40.
14. Rodi W. Comparison of LES and RANS calculations of the flow around bluff bodies // Journal of Wind Engineering and Industrial Aerodynamics 1997. № 6971. P.55-75.
LARGE-EDDY SIMULATION OF FLOW STRUCTURE IN BLUFF-BODY FLAMEHOLDER
© 2013 S.G. Matveev, I.A. Zubrilin
Samara State Aerospace University named after Academician S.P. Korolyov (National Research University)
In this study LES is used in prediction of isothermal flow of triangular bluff-body and compared with prediction provided by URANS simulations and experimental data. Shown, that LES is more accurate in prediction of size of recirculation zone and coherent structures. In addition, lean blowout of flame with using different reaction mechanisms was observed.
Key words: flameholder, computational fluid dynamics, turbulence models, large eddy simulation.
Sergey Matveev, Candidate of Technical Science, Associate Professor at the Heat Engineering and Heat-Engine Department. E-mail: [email protected]
Ivan Zubrilin, Assistant Lecturer at the Heat Engineering and Heat-Engine Department. Email: [email protected]
