CC BY
73
2006 НАУЧНЫЙ ВЕСТНИК МГТУ ГА №108
серия Эксплуатация воздушного транспорта и ремонт авиационной техники. Безопасность полётов
УДК 621.452.3
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТОНАЛЬНОГО ШУМА ВЕНТИЛЯТОРОВ СОВРЕМЕННЫХ АВИАЦИОННЫХ ДВИГАТЕЛЕЙ
М.А. НЮХТИКОВ, В.В. СГАДЛЕВ, А.А. РОССИХИН
Статья представлена доктором технических наук, профессором Елисовым Л.Н.
В статье рассмотрены основы численного метода моделирования акустических характеристик современных авиационных двигателей.
Ввиду высокой важности обеспечения высоких акустических характеристик для авиационных двигателей нового поколения становится актуальным вопрос моделирования их акустических характеристик на раннем этапе проектирования. Реактивный двигатель является сложным источником шума, поскольку шум образуется во всех его узлах: компрессоре, вентиляторе, камере сгорания, турбине, реактивном сопле и струе. Однако доминирующим источником шума для двигателей с большой степенью двухконтурности является тональный шум вентилятора, вызванный ротор-статор взаимодействием и шумом ударных волн.
Моделирование ротор-статор взаимодействия обычно проводилось с использованием схем второго порядка на основе уравнений Эйлера или Навье-Стокса [2]. Для адекватного учета особенностей нестационарного течения в вентиляторе при таком подходе необходимо использовать очень подробные сетки, что приводит к очень большим временам расчета и затратам памяти ЭВМ. Необходимость использования очень подробных сеток, а также низкая точность этих схем при решении задач о распространении волн ставит под сомнение их эффективность. Использование же схем высокого порядка существенно снижает размеры расчетных сеток и существенно повышает скорость и точность расчета. Однако время расчета остается все же очень большим.
В отделении 100 ЦИАМ разработан новый численный метод 3D расчета шума вентилятора с учетом ротор-статор взаимодействия. Метод базируется на прямом численном решении нелинейных и линейных невязких нестационарных трехмерных уравнений для возмущений численными методами вычислительной акустики. Причем, если этого требует решаемая задача, уравнения могут решаться в частотной области методом установления по псевдо времени. Это может быть необходимо, например, при решении задач излучения и распространения шума для отдельных частот.
Расчет генерации и распространения шума проводится в 2 этапа:
1. Расчет 3D стационарных полей течения в воздухозаборнике и в ступени вентилятора. Расчет 3D стационарного поля течения в воздухозаборнике проводится с использованием полной системы уравнений Навье-Стокса и полуэмпирической модели турбулентности в приближении «mixing plane.»
2. Нестационарный 3D расчет распространения возмущений. Для проведения нестационарного расчета шума взаимодействия используются локальные вращающиеся системы координат, связанные с роторами (или вращающиеся с ротором и неподвижная со статором). На их границах для передачи возмущений использованы интерфейсы для “сращивания” вращающихся сеток. Причем передача возмущений организована так, чтобы обеспечивать сквозное прохождение волн через границы и передавать возмущения, вызванные движением среднего поля течения, вызванного относительным движением систем координат.
Расчет излучения может проводиться одновременно с расчетом взаимодействия, но может проводиться и отдельно от него. В первом случае, к расчетной области добавляется геометрия
воздухозаборника. Во втором случае, расчет излучения проводится отдельно, с использованием в качестве источников данных, полученных из нестационарного расчета.
Для расчета шума в дальнем акустическом поле используется подход Фокс-Вильямс-Хокингса [3]. С использованием данных об источниках на поверхности расчетной области, рассчитываются уровни звукового давления в заранее выбранных точках в дальнем акустическом поле.
Дискретизация по пространству построена на основе метода конечных объемов, и с использованием обобщенной на метод конечных объемов DRP (Dispersion Relation Preserving) схемы. Для дискретизации уравнений по времени использована оптимизированная 4-шаговая схема Рунге-Кутта (LDDRK - Low Dissipation and Dispersion Runge-Kutta Schemes) второго порядка, что обеспечивает высокую точность схемы при решении нестационарных задач.
Основная система уравнений
Система линеаризованных уравнений Эйлера в консервативной форме может быть получена из системы полных уравнений Эйлера, записанных в консервативной форме, путем разделения параметров потока на среднее стационарное поле течения (обозначено индексом 0) и нестационарное поле пульсаций (обозначено штрихом):
Р = Ро +Р' Р = Ро + Р '
u = u0 + u
V = v0 + v w = w0 + w ' e = e0 + e '
Подставляя их в систему уравнений Эйлера, записанную в дивергентной форме, и выбрасывая члены второго и более порядка малости, получим линеаризованные уравнения Эйлера:
дО' дЁ' дР' дв
- +---L +----L + -
dt
где
U '
Р'
р0 u '+ u0 • Р' Ро -v'+ V Р' р0 •w '+ w0- Р'
e '
дх
С
: 0
Е1
dy dz
Р0 •u '+ u0 •Р '
(1)
Л
2 • Р0 u • u '+ u0 • Р '+ p '
Р0 • u0 ^'+Р0 u '• V0 +РКu0 •V0
Р00 •u0 •w '+Р0 u ' w0 +РК u0 •w0
( e0 + Р0 )•u '+( e '+ Р ') ^0
для Р \ , в \ - аналогично.
В этих выражениях возмущение энергии е' записывается в следующем виде:
+ (Г +Р')>0 и'+ У0-V'+ ^') +
g-1
2' Р'(ио + V + ™02) + 2' Р'(и '2 + У '2 + ™ ’2)
Для вращающейся с угловой скоростью ю вокруг оси X системы координат система линеаризованных уравнений Эйлера записывается в следующем виде:
dU' дЕ ' д
-------1--------1---
dt дх
F w • zU ' д G '+№• zU '
+
дУ
dz
0 0
где: I = о-и4'
-о-из 0
Здесь из' и и4' определяются в соответствии с выражением для и1' при 1=3 и 1=4.
В том случае, когда исследуемое возмущение может быть охарактеризовано конечным набором частот, нестационарную задачу можно свести к набору стационарных. Так как полученная система уравнений (2) линейна по вектору и', поэтому мы можем разложить вектор и' по времени в ряд Фурье в каждой точке:
и'(х, у, г, і) = и (X, у, 2) + 2 ип (X, у, 2) - е~1'°>п 1
п
где и0 (х, у, г) - 0 гармоника ряда Фурье; ип (х, у, г) - п-я гармоника ряда Фурье. Наша система перепишется в виде:
(3)
дЁ„' д
дх„
f '-a-z-u.
д
дУг
+-
оя'+п-y-un
dzr
■L + i - w -U„
(4)
Для решения полученной системы уравнений будем использовать метод установления по псевдо-времени [4], [6]. Для этого добавим фиктивную производную по времени к нашей системе уравнений:
дЁА+дЁ1+д.
dt дх„
[ Fn '-a-z-Un •] д [о„ '+a-y-Un •]
дУг
dz„
— F + i - w - U„' — S
(5)
Пространственная дискретизация
За основу использованной здесь разностной схемы был взят DRP (Dispersion-Relation Preserving) подход, предложенный К.Тамом [1]. Там и Вебб оптимизируют коэффициенты конечно-разностной схемы, с тем чтобы минимизировать эффекты диссипации и дисперсии. Их подход может быть применен и к методу конечных объемов, как в нашем случае. Так как DRP схемы оптимизируются на равномерной сетке, то нужно преобразовать физическую систему координат (X,Y,Z) в систему координат (I,J,K), где сетка будет равномерной, прямоугольной и иметь размер ячеек равный 1 по всем 3-м осям:
\ х — X (I, J, К ) у — Y (I, J, K ) z — Z (I, J, K )
Тогда система уравнений (1) в новой системе координат запишется в виде:
ду дА дВ дС .
-?- +----+-----+------— 0, (2)
дї ді дJ дК
где
— u1
У det(3),
А, В, С определяются через потоки системы (1) и Якобиан.
Матрица 3 - матрица Якоби, которая вычисляется численно, и имеет вид:
' д1 д1 д1 Л
дх дУ дг
dJ dJ дJ
дх дУ дг
дк дк дк
к дх дУ дг )
Для решения описанной выше системы уравнений используются оптимизированные схемы высокого порядка:
• дискретизация по пространству построена на основе метода конечных объемов и с использованием обобщенной на метод конечных объемов DRP (Dispersion Relation Preserving Scheme) схемы;
• для дискретизации уравнений по времени использована оптимизированная 4-шаговая схема Рунге-Кутта (LDDRK - Low Dissipation and Dispersion Runge-Kutta Schemes) второго порядка.
Результаты расчета
Общий вид вентилятора с внешней обечайкой приведен на рис. 1 и без верхней обечайки на рис. 2. Поле среднего течения рассчитано с использованием программного комплекса ЦИАМ 3DFS. Оно было получено как численное решение усредненных по Рейнольдсу уравнений На-вье-Стокса с моделью турбулентной вязкости nt90. Задача решалась методом Годунова 2-го порядка по пространственным переменным. При этом использовались следующие граничные условия:
На внешней границе расчетной области ставились условия затопленного пространства. На выходе воздухозаборника задавалась величина давления, обеспечивающая заданный расход газа. Расчет производился методом установления с локальным шагом по времени. Всего было сделано 1000 итераций. На рис. 3 приведено поле числа Маха в расчетной области.
Расчет акустических характеристик вентилятора в передней полусфере проводился программным комплексом ЦИАМ 3DAS для основной гармоники с частотой шума ударных волн. Результат был получен как численное решение линеаризованных уравнений Эйлера для осевых гармоник в частотной области. Расчет проводился методом установления по псевдо-времени, для ускорения сходимости был использован локальный шаг по времени (т.е. шаг по времени определялся исходя из условия устойчивости численной схемы для каждой ячейки в отдельности и был разным для всех ячеек). В качестве начального приближения было взято поле, где все величины пульсаций (как действительная часть, так и мнимая) были положены равными 0.
В качестве источника тонального шума на частоте ударных волн использовалось собственное решение для акустических волн в канале с азимутальным числом 15, что соответствует числу лопаток вентилятора, и радиальным числом 1. Амплитуда бралась из расчета ротор-статор взаимодействия.
В результате расчета были получены поля действительной и мнимой части статического давления для всех указанных гармоник. Для того чтобы получить сошедшееся решение потребовалось сделать 20000 итераций. График сходимости нормы решения приведен на рис. 4. Поле пульсаций статического давления в ближнем поле приведено на рис. 5. Видно, что полученная картина излучения похожа на ту, которую мы обычно наблюдаем в эксперименте.
Для расчета излучения воздухозаборника в дальнем поле использовался метод Фокс Вильямса, встроенный в программный комплекс ЦИАМ ЗЭЛБ. При расчете пульсаций давления регистрация давления производилась в точках, выбранных таким образом, чтобы они совпадали с точками регистрации натурного эксперимента.
В процессе расчета использовалась поверхность интегрирования, показанная на рис. 5.
В качестве исходных данных для наших вычислений были взяты результаты расчета акустических характеристик вентилятора в ближнем поле, полученные на первом этапе работы. На рис. 6 проведено сравнение диаграмм направленности для частоты следования лопаток с результатами, полученными по полуэмпирической методике ЦИАМ. Видно, что диаграммы направленности хорошо совпадают в районе максимума и имеют существенное отличие в диапазоне малых углов. Расхождение может быть объяснено тем, что в расчете мы использовали только 1 гармонику, в то время как в действительности из-за наличия неравномерности течения на входе в воздухозаборник и ошибок в угле установки лопаток мы имеем множество гармоник на той же частоте, но с различными азимутальными числами, шум от которых заполняет этот пробел.
Выводы
В работе проведен численный расчет излучения шума вентилятора из воздухозаборника в ближнем и дальнем поле.
• В ближнем поле в результате расчета, основанного на использовании статических уравнений, были получены поля действительной и мнимой части статического давления на частоте следования лопаток.
• Получена диаграмма направленности в дальнем поле, на основе метода Фокс Вильямса. Диаграмма направленности имеет хорошее совпадение с диаграммой направленности, рассчитанной по полуэмпирической методике в районе максимума, и значительное отличие при небольших углах.
Для установления причин расхождения диаграмм направленности на небольших углах необходимы дополнительные исследования. Особое внимание, на наш взгляд, следует обратить на тщательный анализ гармонического состава возмущений на входе в вентилятор, которые вызваны возмущениями среднего поля течения на входе в воздухозаборник.
Рис. 1. Трехмерная математическая модель Рис. 2. Трехмерная математическая модель
вентилятора вентилятора без верхней обечайки
0.1
0.01
Рис. 3. Поле числа Маха среднего течения
0.0001
1e-005 -
Re(<p>)
Im(<p>)
0
4000 8000 12000 16000 2
N
Рис. 4. График сходимости нормы решения при расчете излучения шума в частотной области
80
100
Рис. 5. Поле пульсаций давления для первой распространяющейся гармоники в сечении, проходящем через центр расчетной области и ось воздухозаборника
Рис. 6. Диаграмма направленности суммарного уровня звукового давления
ЛИТЕРАТУРА
1. Tam C. K. W., Webb J. C. Dispersion -Relation-Preserving Schemes for Computational Acoustics, // Journal of Computational Physics, 107, 262-281, 1993.
2. Hu F. Q., Hussaini M. Y., Manthey J. Low-Dissipation and -Dispersion Runge-Kutta Schemes for Computational Acoustics. NASA Report, 1994.
3. Williams J.E., Ffowcs Hawkings D.L. Sound Generated by Turbulence and Surfaces in Arbitrary Moution // Philosophical Transactions of the Royal Society. Vol. A264. No 1151, 1969. PP. 321-342.
NUMERICAL SIMULATION OF ACOUSTICAL CHARACTERISTICS OF MODERN JET
ENGINES
Nyukhtikov M.A., Sgadlev V.V., Rossikhin A.A.
This paper present the basis of the numerical method for modeling acoustical features of modem engines used in
CIAM.
Сведения об авторах
Нюхтиков Михаил Анатольевич, 1976 г.р., окончил Московский физико-технический институт (1999), кандидат физико-математических наук, начальник сектора вычислительной акустики лопаточных машин ЦИАМ, автор более 10 научных работ, область научных интересов - численное моделирование акустических характеристик.
Сгадлев Василий Викторович, 1983 г.р., окончил МГТУ ГА (2005), инженер ЦИАМ, область научных интересов - численное моделирование акустических характеристик.
Россихин Антон Анатольевич, 1977г. р., окончил МГУ (2000), кандидат физико-математических наук, младший научный сотрудник ЦИАМ, автор более 10 научных работ, область научных интересов -вычислительная акустика.
