CC BY
87
Том XL V
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
2014
№ 2
УДК 534.231 519.63 629.7.035
РАЗРАБОТКА МЕТОДОЛОГИИ РАСЧЕТА ШУМА ВИНТОВ С ИСПОЛЬЗОВАНИЕМ СУПЕРКОМПЬЮТЕРОВ
В. Ф. КОПЬЕВ, В. А. ТИТАРЕВ, И. В. БЕЛЯЕВ
Описывается методология расчета аэродинамических и акустических характеристик вращающегося винта реальной формы. Численный метод расчета аэродинамических характеристик основан на использовании неявной схемы высокого порядка аппроксимации на произвольных неструктурированных сетках. Звуковое поле вычисляется с использованием ФВХ-подхода. Для ускорения счета реализовано распараллеливание решателя на основе технологии MPI. Работоспособность метода и комплекс программ иллюстрируются на задаче обтекания шестилопастного авиационного винта и включают валидацию решателя путем сравнения данных расчета с экспериментом, тестирование масштабируемости метода до 1024 процессорных ядер, исследование влияния величины шага по времени на скорость сходимости к стационарному решению. Расчет акустического поля верифицирован на основе точного решения линейного приближения.
Ключевые слова: винт, неструктурированная сетка, ТВД-схема, суперкомпьютерные расчеты, шум обтекания.
ПРИНЯТЫЕ ОБОЗНАЧЕНИЯ
Аа — 1-я грань г'-й ячейки сетки; С — число Куранта;
с0 — скорость звука в невозмущенном потоке;
gy, — компоненты метрического тензора;
Г, С, Н — конвективные потоки в декартовой системе координат; Н — толщина лопасти винта; кп — пП0/с0, волновое число; М0 — и о/ с0, число Маха;
КОПЬЕВ Виктор Феликсович
доктор физико-математических наук, начальник отделения ЦАГИ
ТИТАРЕВ Владимир Александрович
кандидат физико-математических наук, ведущий научный сотрудник ЦАГИ и ВЦ им. А. А. Дородницына РАН
БЕЛЯЕВ Иван Валентинович
кандидат физико-математических наук, старший научный сотрудник ЦАГИ
Mл — и л/ с0, число Маха; п — номер Фурье-гармоники;
п — внешняя нормаль к возмущенной поверхности лопасти винта; п0 — внешняя нормаль к невозмущенной поверхности; — общее количество ячеек сетки;
р — давление возмущенного потока; р0 — давление невозмущенного потока;
г, ф, 0 — сферическая система координат; Я — радиус винта; К — вектор невязки; Г — время; и — Я/¥„, характерный масштаб времени; и, V, — компоненты скорости течения, создаваемые винтом в точке поверхности ФВХ, в системе координат наблюдателя; и — вектор консервативных переменных; и0 — скорость поступательного движения винта;
vb v2, v3 — компоненты скорости движения точки поверхности ФВХ в системе координат наблюдателя;
Vx, Vy, Vz — компоненты скорости вращения системы координат;
V* — с0/^/то", характерный масштаб скорости;
Vol; — г'-я ячейка сетки;
|Vol;| — объем г'-й ячейки;
W — (ра , u, v, w, p) , вектор простых переменных;
x1, x2, x3 — координаты точки положения наблюдателя в системе координат наблюдателя;
x, y, z — декартова система координат;
х, y, z — движущаяся система отсчета, связанная с осью винта;
У1, У2, У3 — координаты точки положения источника в системе координат наблюдателя;
Y о — отношение удельных теплоемкостей;
П1, П2, П3 — пространственная координата точки поверхности ФВХ в системе отсчета источника;
р, Е, у — винтовая система координат;
р, ф, z — цилиндрическая система координат;
р0 — плотность невозмущенного потока;
р а — плотность возмущенного потока;
т — временная координата точки поверхности ФВХ в системе отсчета источника;
ф — потенциал возмущенного течения;
юп — пПо, частота;
ю* — 1/1*, характерный масштаб частоты;
Q0 — угловая скорость вращения
U Л
1. ВВЕДЕНИЕ
В настоящее время различные научные коллективы в основных авиапроизводящих странах, в том числе в России, активно включились в разработку новых схем авиационных двигателей. В качестве альтернативы современным ТРДД для обеспечения большей экономичности, в частности, рассматриваются ТВВД с открытым винтовентилятором различных схем (одновинтовая схема, схема с соосными тянущими или толкающими винтами). Компоновка с открытым ротором рассматривается в качестве серьезной перспективы для силовых установок узкофюзеляжных самолетов, которые придут на смену самолетам семейства Airbus 320 и Boeing 737. Учитывая исключительную актуальность проблемы шума винтов при создании этой технологии, в акустиче-
ском отделении ЦАГИ разрабатывается комплекс программ численного моделирования задач аэродинамики и акустики винтов, позволяющий исследовать различные конфигурации с точки зрения как шума на местности, так и шума в салоне.
Первые попытки расчета шума винтов относятся к работе Л. Я. Гутина [1], использующей для расчета шума нагрузки интегральные аэродинамические и суммарные геометрические характеристики винта. Учет шума вытеснения, связанного с толщиной лопасти [2, 3], позволил более точно рассчитывать акустическое поле винтов, и после существенной доработки [4, 5] такой подход лег в основу многих полуэмпирических методик [6, 7]. Дальнейшее развитие аналитических подходов позволило учитывать распределенные дипольные источники, связанные с локальными силами воздействия лопастей на поток; монопольные источники, связанные с вытеснением движущимися лопастями объема газа или жидкости, а также квадрупольные источники (в основном связанные с нелинейными возмущениями потока в окрестности высокоскоростных лопастей) [8—10]. Значительный прогресс был достигнут в работах Хансона [11 —13], использующих для расчета звукового поля уравнение Фокс Вильямса — Хоукингса [14].
В последующие годы методы расчета шума винтов получили существенное развитие, связанное в первую очередь с использованием новых вычислительных мощностей, имеющихся в наличии у научных коллективов [15]. Современные подходы [16—18] позволяют вычислять не только аэродинамические и акустические характеристики современных винтов, но и влияние на них эффекта установки винта на самолете, неоднородности внешнего потока и т. д.
В настоящей работе на примере расчета обтекания шестилопастного авиационного винта (рис. 1) рассмотрен комплекс вопросов, направленных на создание собственного вычислительного алгоритма для расчета нагрузок, тяги и акустического поля винта самолета. Отличительными особенностями создаваемого комплекса программ являются поддержка форматов неструктурированных сеток ведущих коммерческих сеткопостроителей (Neutral пакета Gambit и StarCD пакетов ICEM CFD и Gridgen), реализация явной и неявной схем высокого порядка аппроксимации на произвольных неструктурированных сетках, эффективное распараллеливание решателя на основе технологии MPI для использования на современных суперкомпьютерных системах, а также вывод данных, адаптированных для акустических программ ЦАГИ и ведущих пакетов визуализации. Использование разрабатываемого вычислительного комплекса программ позволяет создать среду моделирования, включающую в себя все этапы решения прикладных задач: получение геометрической модели объекта, построение расчетной сетки в физическом пространстве, расчет обтекания с использованием суперкомпьютеров, обработку результатов в акустическом постпроцессоре и выдачу рекомендаций разработчикам.
Рис. 1. Геометрия винта с поворотным механизмом
Расчет аэродинамических нагрузок на лопасти винта дает возможность провести валидацию метода путем сравнения данных расчета с экспериментально полученными интегральными характеристиками тяги и момента. При отсутствии надежных экспериментальных данных по акустическим характеристикам рассматриваемого в работе винта проблема валидации численного метода расчета шума оказывается трудноразрешимой. Однако в случае малых углов отклонения лопастей для математического описания задачи можно использовать линейное приближение, которое дает точные выражения для двух компонент шума винта: шума вытеснения (связанного с геометрией лопасти) и шума нагрузки (связанного с силами, действующими со стороны лопасти на газ). Имея эти данные, результаты численного метода можно сравнить непосредственно с точными выражениями математической задачи линейного приближения и тем самым обеспечить верификацию расчета.
Численный расчет шума в настоящей работе осуществляется путем использования интегрального подхода Фокс Вильямса — Хоукингса (ФВХ), позволяющего пересчитать аэродинамические нагрузки на контрольной поверхности в звуковое поле. ФВХ-подход представляет собой основной инструмент, с помощью которого в современных вычислительных пакетах определяется шум в дальнем поле. В своей общей формулировке данный подход позволяет учесть все источники аэродинамического шума. В качестве поверхности ФВХ, на которой задаются давление и скорость течения, в данной работе использовалась поверхность лопастей винта. В дальнейшем планируется рассмотреть семейство вложенных поверхностей ФВХ, учитывающих распределенные квадрупольные источники в области между лопастью и поверхностью ФВХ, с контролем сходимости по различным поверхностям в акустической зоне. Верификация вычислительного комплекса осуществляется с помощью точного выражения для шума нагрузки и шума вытеснения, полученного в линейном приближении, когда квадрупольные источники заведомо малы.
Таким образом, в работе проведен полный аэродинамический расчет обтекания шестилопа-стного авиационного винта (см. рис. 1) для различных значений скорости вращения и скорости набегающего потока. Все численные результаты получены из решения уравнений Эйлера для сжимаемого газа во вращающейся системе координат и включают валидацию решателя путем сравнения интегральных данных расчета с экспериментом, тестирование масштабируемости метода до 1024 процессорных ядер, исследование влияния величины шага по времени на скорость сходимости к стационарному решению, верификацию численного решения по точному решению задачи линейного приближения.
2. ИЗЛУЧЕНИЕ ШУМА ВИНТОМ В ЛИНЕЙНОМ ПРИБЛИЖЕНИИ
Вычисление шума винта в линейном приближении направлено на проведение верификации решения, полученного развитым в работе численным методом, позволяющим считать аэродинамические и акустические характеристики винтов с высокой точностью. Приведенное в настоящей работе аналитическое решение [19], в отличие от принятого в литературе подхода Хансона, использующего интегральное уравнение Фокс Вильямса — Хоукингса, получено прямым использованием граничных условий на поверхности лопасти, что делает вывод основных соотношений наиболее наглядным. Это существенно в контексте настоящей работы, поскольку задача о шуме винта в линейном приближении применяется здесь для верификации численного метода.
В настоящем разделе определены линейные составляющие шума винта — шум вытеснения и шум нагрузки, представленные в виде рядов Фурье в движущейся системе координат. Выражения справедливы как в ближнем, так и в дальнем поле винта, при любой скорости обтекания. Определенная громоздкость математического аппарата, используемая при выводе основных соотношений, не должна закрывать элементарную сущность вывода: решаются линеаризованные уравнения сжимаемого газа с условием непротекания на поверхности лопасти. Эти условия непосредственно приводят к нахождению плотности распределенных источников в случае шума вытеснения и напрямую связаны со скачком подъемной силы в случае нагруженной лопасти. В линейном приближении эти два источника очевидно независимы.
2.1. ВИНТОВАЯ СИСТЕМА КООРДИНАТ
Рис. 2. Цилиндрическая система координат
Свяжем исходную систему отсчета с землей, т. е. будем считать, что жидкость на бесконечности покоится. Введем цилиндрическую систему координат р, ф, г, ориентированную так, чтобы при своем движении лопасть совершала поступательное движение со скоростью -и0) вдоль г и вращательное с угловой скоростью 00 вдоль ф (рис. 2, а):
x = р cos ф, У = -р sin ф, z = z.
(2.1)
Введем на поверхности каждого цилиндра р = const вместо координат ф и z координаты Е и у по формулам
. U0 р^0 z = y cos a-csrn a = y——-ц--,
U U
U Л U Л
С р^0 , U0 РФ = -Y sin a -С cos a = -y--ц-
U Л U Л
(2.2)
где cos a
= U0/UЛ, sina = pQ0/UЛ, UЛ =yjU02 + р2Ц^ (рис. 2, б). Тогда для С, Y и
р получим:
U0 р^0 ц = -ТГ~ РФ-' 0
U
Л
U
Л
U0 р2Ц,
Y = —— z--- ф,
Uh Uл
Р = Р.
(2.3)
Система координат Е, у, р (2.3) называется винтовой и построена таким образом, чтобы каждая точка лопасти описывала при своем движении геликоидальную линию Е = const,
р = const, -да < y < Система Е, у, р имеет период, поэтому для однозначности на каждой цилиндрической поверхности р = const выбираем, что для лопасти без отгиба Е меняется в пределах (0, 2прио/UA), -да<у<да. Несмотря на кажущуюся простоту, эта система координат имеет довольно сложную пространственную структуру, связанную с зависимостью UA, а, следовательно, и 9, от р. Удобство выбранной системы координат заключается в том, что любая тонкая поверхность без угла атаки, лежащая на поверхности Е = 0 (более общее условие Е = f (р), соответствующее лопасти с отгибом), не создает при своем геликоидальном движении никакого возмущения потока. В этом смысле система координат Е, у, р так же естественна для изучения винтов, как и простая декартова система координат для тонкого крыла.
Для правильной записи векторных соотношений, определяющих граничные условия на лопасти, необходимо иметь выражение для метрического тензора g'J в системе у, р. Вычисляя
значения производных «старых» координат р, ф, г по «новым» у, р и используя выражение для метрического тензора в цилиндрической системе координат, можно получить, что в винтовой системе
' 1 0 - 4 л 1+ АгВг 4"
g.J = 0 1 - В ? ^ = 4В* 1+В2 Вг ?
Г 1" + -ч, 4 V В* 1
= а^Ц =1, 4 = и 2 и0 ТГ 2 и Л р ' и лV р у + и л 2и и л Л Е у
Громоздкость метрического тензора в винтовой системе связана с упоминавшейся выше зависимостью и л от р, приводящей к кручению репера у при изменении р.
2.2. ЛОПАСТЬ БЕЗ УГЛА АТАКИ. ШУМ ВЫТЕСНЕНИЯ
Для лопасти без угла атаки граничные условия на верхней и нижней поверхностях достаточно просто могут быть связаны со скоростью набегающего потока и локальной толщиной лопасти, т. е. чисто математически движение такой лопасти может быть заменено с помощью граничных условий движением распределенных стационарных источников с известной плотностью.
Будем рассматривать лопасть без отгиба со средней линией Е = 0 (случай ¡ = / (р) рассмотрен в Приложении 1) под нулевым углом атаки. Лопасть за счет толщины выступает за пределы невозмущенной поверхности ¡ = 0, поэтому нормаль не совпадает в точности с вектором Е и возмущает поток. Уравнение верхней и нижней поверхности:
¡Г =8/1 (у, р), Г =8/2 ( р)
(2.5)
/ - /2 = Ь (у, р) — толщина лопасти. При рассмотрении шума вытеснения полагаем / =-/2.
Считаем е ^ 1.
Граничное условие непротекания имеет обычный вид:
(-ил +8УФ, я]
= 0,
(2.6)
которое в силу линейности задачи может быть снесено на невозмущенную поверхность Е = 0. Здесь (—ил) — скорость набегающего основного потока; Ф — потенциал возмущенного течения; п — внешняя нормаль к возмущенной поверхности:
я+ =
5у ф
(
п =
-1, е /, е /
5у ф
Л
(2.7)
И1 = g'JnInJ
¡=0
= 1.
Формула (2.6), примененная к верхней и нижней поверхности лопасти, дает выражение для скачка нормальной производной потенциала:
[МФ)>
дФ
дп0
= и .
Е = 0
дк ду.
(2.8)
Здесь п00 — внешняя нормаль к невозмущенной поверхности п00 =(1, 0, 0), и'А =(0, - ил, 0).
Всюду вне этого участка потенциал Ф вместе со своими первыми производными непрерывен и удовлетворяет волновому уравнению. Таким образом, задача о шуме вытеснения сводится к решению обыкновенного волнового уравнения для потенциала с заданным скачком его нормальной производной на заданной поверхности. Для волнового уравнения, как будет показано ниже, скачок нормальной производной [дФ/дп] математически эквивалентен распределению
монопольных источников в правой части. Таким образом, локальность граничных условий позволяет свести задачу о шуме вытеснения к решению обыкновенного волнового уравнения с известной правой частью.
Однако волновому уравнению во всем пространстве удовлетворяет не только потенциал Ф, но и давление р, использование которого оказывается во многих отношениях удобнее. Действительно, для шума нагрузки, который будет рассмотрен ниже (т. е. лопасть под углом атаки), потенциал испытывает скачок не только на самой лопасти, но и на всей винтовой полубесконечной поверхности вниз по потоку, что приводит к необходимости, если задача решается для потенциала Ф, распределять источники на всей полубесконечной поверхности, включая лопасть и след. Этого, очевидно, не происходит для давления р, непрерывного всюду, кроме самой лопасти. Исходя из такого очевидного преимущества для шума нагрузки, целесообразно в случае шума вытеснения перейти от рассмотрения потенциала к рассмотрению давления. Для этого необходимо перейти от граничного условия (2.8) для скачка потенциала [дФ/дп] к граничному условию
на скачок [дФ/дп], что можно сделать, воспользовавшись уравнением Эйлера:
д р дФ р _УФ = -У^_ или -= -—
д р0 д р0
(2.9)
где Р0 — плотность газа на бесконечности.
В системе координат Е, У, Р зависимость от времени, очевидно, входит только в комбина-
д
ции (у + и). Поэтому д\д1 = ил — и (2.9) можно переписать в виде:
ду
и дФ
Р = -Р0ил^".
ду
(2.10)
. | д д д ) Учитывая, что и'л =(0, - ил, 0), Уг- =1 —, —, — , соотношение (2.10) можно записать
^д^ ду др
в инвариантном виде:
р = Р0 (и лУФ).
др
(2.11)
Чтобы получить нормальную производную -= (п00Ур), подействуем на обе части равен-
дп
ства
(2.11) оператором (п0, У):
(п0Ур ) = |р = Р0 (п>У)(и лУФ).
(2.12)
Преобразуем правую часть равенства (2.12). Имеем:
(n0V)(UЛУФ) = (UaV)( п0УФ) + q,
где q =
дФ
dxj
U - u
dx,
dn0 dx,
дФ
dy
(
dU л dp
- U,
dn0 dY
Покажем, что q = 0, т. е. оператор (^У) можно проносить через (ТлУ). Действительно,
с учетом конкретного вида нормали выражение в круглых скобках равно нулю, и (2.12) можно записать в виде:
Ф =Ро (Uл^^Ф^и лдГдФ'
dn,
ду
Vdno J
(2.13)
Используя формулу (2.8), получим известное выражение
дРо
dn0
fl2 d2h ti ч
= -poUл—=q (у, p),
5=0
dY 2
(2.14)
полученное Хансоном в [12]. Для лопасти с отгибом Ф 0 связь (2.14) не имеет места, так как в этом случае нормаль к лопасти уже не будет лежать на поверхности цилиндра p = const. Отметим, что граничное условие (2.14) записано в начальный момент времени. С учетом движения лопасти надо всюду заменить зависимость от у на у + Uлt.
Таким образом, задача о шуме вытеснения сводится к краевой задаче для волнового уравнения
{АР}-1 ^ = 0,
dp dn0
c02 dt2
p u2 d2h = —P0UЛ
(2.15)
^ = 0
ay
2
с условием обращения возмущений в нуль на бесконечности.
Заданному скачку [ф/дщ ] на поверхности лопасти легко сопоставить распределение источников. Воспользуемся известной формулой [20, стр. 119]
А/ = {А/} + £ A ([/ ] cos nXI 5s) + £ [{{/dx, } ] cos nXI 8S,
i=1 dXi i=1
(2.16)
где {А/} означает обычный лапласиан функции вне линий разрывов; 55 — простой слой на поверхности разрывов 5; [/] и ф } ^ — скачки функции и ее производных через поверхность 5. Формула (2.16) понимается в обобщенном смысле. В нашем случае формула (2.15) перепишется в виде:
{Ар} = Ар -
dp dn
5s.
(2.17)
Тогда, заменяя в волновом уравнении (2.15) {Ар} согласно выражению (2.17), получим вместо краевой задачи обобщенную, с известными источниками в правой части:
Ар -
_L dip
c0 dt2
dp dn0
5s - q'(y + илt, p)5(( = -p0 иЛ ^ 5(0.
(2.18)
Формула (2.18) является главной формулой для задачи о движении лопасти без угла атаки в сжимаемом газе. При ее выводе использовалось только предположение о малости возмущений, во всем остальном она является точной и справедлива при любых относительных скоростях движения лопасти.
Решение уравнения (2.18) представляет собой свертку правой части с функцией Грина волнового уравнения [3]. Однако для получения окончательного результата удобнее перейти в движущуюся систему отсчета х, у, 5, связанную с осью винта:
у = у,
(2.19)
г = г + U0t.
Именно в этой системе отсчета, движущейся со скоростью —и0 вдоль оси г, источник становится периодической во времени функцией, и для получения решения можно пользоваться разложениями в ряд Фурье. Обычное волновое уравнение при этом перейдет в конвективное. С учетом формул (2.3) и (2.19), источник перепишется в виде:
г(ф—М, рУц1 р(ф—М+и
и Л
и Л
У V ил
(2.20)
л У
Так как процесс периодический с периодом Т = 2п/О0, можно решать задачу для каждой Фурье-гармоники отдельно. Уравнение (2.18) имеет вид:
А —I -кп +Мо —
Рп = Чп, кп =
п О0
Мо =Ц0,
(2.21)
где
)е'п Оот а т,
О 2п/"о
Чп (^ ^ = 1л ^ Ч *, т)
о
О 2ПО0
Рп (р, ф) = О° Г р(р, Ф, т) е'пОотат,
(2.22)
коэффициенты разложения ч и р в ряд Фурье
да да
Ч = £ чпе~тО°, р = £ рпе~1пО°. (2.23)
Обратим внимание, что для каждой гармоники источник звука является объемным распределением, с постоянной амплитудой в каждой точке объема, даваемой выражением (2.22), и занимает всю область (р, 5), ометаемую винтом, а не расположен на лопасти и, тем более,
не вращается вместе с ней. Источник звука есть правая часть в волновом уравнении (2.21) и имеет ту структуру, которая получается после сведения задачи к волновому уравнению с правой частью. Вращающаяся сила или вращающийся объем, эквивалентные вращающимся аэродинамическим источникам (2.8) или (2.14), есть только причины того, что в задаче появляются источники звука, структура которых дается (2.22).
Решение уравнения (2.21) для каждой отдельной гармоники не зависит от того, какой режим обтекания возникает на концах лопасти — дозвуковой или сверхзвуковой. Наличие транс-
звуковых зон, а также ударных волн, возникающих при сверхзвуковом обтекании, конечно, ограничивает применение линейной теории. Однако мы исходим из того, что для относительно тонких лопастей, как правило, применяемых в случае больших скоростей, акустическое излучение определяется в основном только линейными эффектами. В этом случае необходимо выбрать просто причинную функцию Грина уравнения (2.21), соответствующую уходящим на бесконечности волнам. Сложности в случае сверхзвукового режима возникнут при суммировании рядов (2.23). Очевидно, что наличие слабых разрывов (местных скачков) приведет к медленному убыванию коэффициентов ряда (2.23) и необходимости суммирования большого числа членов.
Для источника звука вида (2.20) легко вычислить в явном виде коэффициенты дп:
и Л
2пЦ0р
ц
и о
р
ио
(2.24)
Эта запись означает, что в формуле (2.14) для Ч надо у заменить на
Цл ио
-г.
2.3. ЛОПАСТЬ ПОД УГЛОМ АТАКИ. ШУМ НАГРУЗКИ
С шумом нагрузки можно поступить аналогично шуму вытеснения. Наличие подъемной силы означает, что на поверхности лопасти в каждой точке имеется скачок давления Ро (у, р) (чтобы тяга была направлена вперед, Р0 должно быть меньше нуля). Эта величина, однако, не может быть просто связана с локальной скоростью набегающего потока и углом атаки, как в случае шума вытеснения, и, строго говоря, вообще не может быть определена отдельно от создаваемого лопастью звукового поля. Для правильного определения нагрузки Ро (у, р) необходимо
решение сложного интегрального сингулярного уравнения.
Однако, поскольку в п. 3. проводится полный аэродинамический расчет поля винта, будем считать нагрузку Ро (у, р) заданной. Тогда по формуле (2.16) такому скачку давления на лопасти сопоставим в уравнении (2.18) плотность вращающихся с лопастью диполей с интенсивностью
ч" = дТС55) = Ро(у, р)5'0
(2.25)
В системе, связанной с осью винта, эта функция будет периодической. Поступая аналогично шуму вытеснения, получим для каждой Фурье-гармоники:
]_ 2п
иоР.
А —I — ¡кп +Мо —
гпРо
Цл V ио
Рп =
^ Р
Р Оо Фо
и л ду
Цл
V ио
^ Р
■ (Оо -гп I —- г
, V ио
(2.26)
где источник звука ч'п локализован в области, занятой лопастью.
2.4. ИЗЛУЧЕНИЕ ЗВУКА ВИНТОМ В СВОБОДНОМ ПРОСТРАНСТВЕ
Таким образом, движущийся винт представляет собой распределенный по объему источник звука, интенсивность которого была определена в предыдущем разделе. В системе координат, движущейся вместе с винтом, этот источник является периодическим, и, следовательно, звуковое поле от него может быть представлено в виде суммы гармоник рп (2.23). Таким образом, задача сводится к определению звукового поля от стационарного распределения гармонических источ-
ников, распределенных по объему, заметаемому лопастями при их вращении. Будем использовать сферические координаты г, ф, 9, цилиндрические координаты р, ф, г и декартовы координаты х, у, г. Оси обеих систем координат направлены параллельно потоку. В движущейся системе координат звуковое поле гармонического источника описывается конвективным волновым уравнением (2.21) или (2.26). Решение этого уравнения представляет собой свертку источника q с функцией Грина О:
Рп (?) = {Оп (г - г") qn (Г)а?,
(2.27)
О( - г') =
„~1кпЯ
4пЯ
(2.28)
ю
кп , Я =
М0 (г'- г) + Я* в2
= л 1 - М,
0
Я* = \/(г'- г)2 +Р2 (X'- X)2 +(у'-у)
Суммарное звуковое поле от источника будет однозначно определено в случае удовлетворения условию излучения при г ^<х>. При представлении суммарного поля в виде ряда по гармоникам условие излучения обеспечивается удовлетворением этого условия для каждой гармоники, что соответствует выбору знака кп/юп > 0 в функции Грина (2.28).
При дозвуковой поступательной скорости винта и о каждая из гармоник представляет собой гладкую функцию, так как определяется решением уравнения (2.21) или (2.26) эллиптического типа. В случае, когда лопасти имеют дозвуковую концевую скорость, суммарное акустическое поле также является гладким, и, следовательно, особенности излучения винта определяются небольшим числом гармоник. В случае, когда лопасти имеют сверхзвуковую концевую скорость, суммарное акустическое поле имеет разрывный характер, несмотря на то, что каждая гармоника непрерывна. Поэтому представление такого поля в виде ряда по гладким гармоникам менее эффективно в связи с медленной сходимостью ряда. Вместе с тем, если требуется найти не временное представление акустического поля винта, а его спектральный состав, то рассматриваемая процедура одинаково эффективна в обоих случаях.
Представление звукового поля в виде свертки (2.27) позволяет вычислить звуковое поле на любом расстоянии от источника, как в дальней, так и в ближней области. В дальней области можно получить асимптотические выражения для звукового поля, раскладывая функцию Грина О(г - г') в интеграле (2.27) по малому параметру в = г'/г. В частности,
Я* = г^ 1 - М2 эш2 9 - 2вЛ* + в2В*
= г
Ф - М2 эт2 9 -
вЛ
Ф - М2 эш2 9 2
- Л
*2
1 - М2 эт2 9
В
) 4
1 - МО эШ2 9
О
(в3)
где
Л* = 008 9 008 9' + Р ЭШ 9 ЭШ 9' 00Э (ф-ф'),
В* = 1 - М2 эШ2 9'.
0
Учитывая, что для точек источника 0 ' ~ п/2, получим:
A ~ р2 sin 9 cos (ф-ф' ),
в* ~i-м2 =р2.
Тогда, оставляя только линейные члены разложения по в, получим для функции Грина (2.28):
-ibknp ' cos (ф-ф')+ia—1 knz'
Gn ( - Г) = G()е Mo , (2.29)
G (r ) = -
3^/l-í
ik„ra\ 11-M2 sin2 f
4nr^1 - M2 sin2 0
a = -
1
1 - M,
M0cos 0 - M2 sin2 0
b =
sin0
■y/i - MO sin2 0
Для источника с плотностью вида дп (г ) = д(р, г)егпф, где п — целое число, подставим (2.29) в (2.27) и вычислим интеграл по ф'. В результате получим в области дальнего поля:
ia—1 kz'
Pn (?) = G(i?)ine'^jj2пе Mo q(p ', z ') Jn (bknp ' )p dp dz',
(2.30)
где Jn — функция Бесселя. Это выражение получено с учетом отбрасывания первого члена в А*,
что может ухудшить точность при малых 0.
В частности, для источников шума вытеснения и нагрузки, плотность которых определяется из (2.24), (2.26), получим п-ю гармонику дальнего поля
п
тф-т—
Рп (Г) = Рп (0)О(г )е ■ 2, где для шума вытеснения и шума нагрузки получим соответственно:
ЕТп (0) = |а%(аЫЛкп)) (Ькпр' )р ', \2 '
(2.31)
fj
(9 ) = -[ in-
U л
иop'
, (p '»о У
1 -a---L-
U Л
po (a M л kn )Jn (bknP ')dP\ мл= U л/ V
Формула (2.31) для дальнего поля удобна при анализе влияния конструктивных параметров винта на уровень его шума. Вместе с тем, вычисление звукового поля для винта заданной конструкции целесообразно проводить на основе точной формулы (2.28), пригодной на любых расстояниях от источника.
Выражение (2.27) с источниками (2.26) и (2.24), полученное в линейном приближении для малонагруженных винтов, представляет собой главное выражение для верификации численного метода расчета шума винтов.
3. ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ЗАДАЧИ ОПРЕДЕЛЕНИЯ АЭРОДИНАМИЧЕСКИХ И АКУСТИЧЕСКИХ ХАРАКТЕРИСТИК ВИНТА
3.1. РАСЧЕТНЫЕ УРАВНЕНИЯ
Рассмотрим обтекание вращающегося шестилопастного винта равномерным потоком воздуха с заданными давлением р0, температурой Т0 и скоростью и0. Скорость вращения винта постоянна и равна N об/с. Геометрия винта представлена на рис. 1 и включает в себя бобышку, поворотные механизмы и лопасти. В некоторых расчетах также использовалась упрощенная геометрия без поворотного механизма. Выберем в качестве масштабов давления, температуры и плотности значения на бесконечности, в качестве масштаба длины — радиус винта Я. В качестве масштабов скорости, времени и угловой скорости примем
V* = С0/, t* = R/V*, ю* = 1/t*,
здесь у0 — показатель адиабаты (отношение удельных теплоемкостей). Далее безразмерные переменные будут обозначаться теми же буквами, что и размерные.
В векторном виде сжимаемые уравнения Эйлера в системе координат, вращающейся вокруг оси z с угловой скоростью Qo, имеют вид [21, 22]:
Uf + Fx + G у + Hz = S, (3.1)
где вектор консервативных переменных U, конвективные потоки F, G, H и источниковый член S задаются выражениями:
U = (р fl, Р fiu-> Р flv, Р flw, E) , S = -Р fl ( - v Q0, + u ^ 0, 0) ,
(3.2)
F = F - Vx U, H = H - Vy U, G = G - Vz U.
Здесь F, G, H — обычные эйлеровские конвективные потоки, вектор (Vx, Vy, Vz) характеризует вращение системы координат с угловой скоростью Qo:
Vx = -y Qo, Vy = x Qo, Vz = 0, (3.3)
pfl — плотность газа, (u, v, w) — поле скоростей газа в системе координат (x, y, z).
Уравнения (3.1), (3.2) используют компоненты скоростей газа в неподвижной системе координат. Решение задачи также можно строить, используя скорость газа относительно вращающейся системы координат [22]. При этом векторы конвективных потоков имеют тот же вид, что и для уравнений в неподвижной системе координат; это позволяет без внесения изменений использовать все имеющиеся численные методы дискретизации конвективных потоков. Однако источниковый член в такой формулировке включает в себя силу Кориолиса и центробежную силу, линейно пропорциональную радиальной координате. Наличие такого источника требует построения сбалансированной (так называемой well-balanced) схемы, которая бы могла сохранять равномерный набегающий поток с радиальной компонентой относительной скорости, растущей линейно пропорционально радиусу. В [22] отмечается в целом меньшая надежность методов решения, использующих относительные скорости, по сравнению с подходом (3.1), (3.2), основанным на использовании абсолютных скоростей.
3.2. ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ЗАДАЧИ 3.2.1. Общая структура метода
Стационарное решение задачи строится методом установления по времени. Введем в физические переменные x = (x1, x2, x3 ) = (x, y, z) расчетную сетку из ячеек Vol;-. Каждая из ячеек может быть тетраэдральной, пирамидальной, гексаэдральной или призматической формы и образована несколькими треугольными или четырехугольными гранями Лц. Общее количество ячеек
равняется Ntot. Объем ячейки обозначим через |Vol;-1. Пусть U" — среднее значение вектора
консервативных переменных в ячейке с индексом i в момент времени tn. Неявная конечно-объемная схема для решения системы уравнений (3.1) записывается в виде:
AU
:R,"+1, AU,. = Uf+1 - Uf, R "+1
At | Vol, | l
Vr Z ф"1 + Sf. (3.4)
Здесь численный поток через грань I ячейки / дается формулой
Фг7 = + пуО + пг н) М, (3.5)
где (пх, пу, п2) — внешняя единичная нормаль к грани. Проводя линеаризацию по времени и записывая разностную схему для всех ячеек сетки, получаем:
Гт Л dRn Л
I -At-
dU
AU = AtRn. (3.6)
Описание численного метода будет законченным, если будут указаны следующие части численного алгоритма: способ вычисления правой части схемы (3.6) с высоким порядком аппроксимации и метод решения нелинейной системы уравнений для приращения решения Аи 1.
3.2.2. Процедура реконструкции
Нахождение численных потоков требует определения средних значений компонент вектора и на гранях ячейки по известным средним значениям в ячейках пространственной сетки. Для получения приемлемой точности счета значения на гранях должны аппроксимироваться со вторым порядком точности и выше. Используемый в настоящей работе алгоритм расчета значений и на гранях ячейки с помощью пространственной реконструкции основан на [23—26]. Ниже приводится его краткое описание.
Для построения схемы используется кусочно-полиномиальное представление каждой компоненты вектора простых переменных W = (р^, и, V, w, р) в пространственной ячейке, обеспечивающее второй порядок аппроксимации по пространству. Ниже для простоты будем выписывать расчетные формулы для скалярной величины /; в вычислениях полученные выражения применяются для каждой компоненты W. Зависимость от времени ^ опущена, поскольку процедура реконструкции выполняется в фиксированный момент времени, соответствующий шагу метода интегрирования по времени.
В каждой ячейке Уо1г- решение локально представляется в виде реконструкционного многочлена рI (х), представляющего собой разложение по базисным функциям. Коэффициенты данного разложения вычисляются при помощи использования значений / в ячейках шаблона реконструкции, который строится путем рекурсивного добавления к ячейке Уо1г- ее соседей до тех пор, пока не достигнуто требуемое число ячеек М в шаблоне. При этом требуется, чтобы в каждой ячейке Уо1т шаблона среднее значение р1 равнялось среднему значению /ш в этой ячейке. Получающаяся переопределенная система уравнений решается методом наименьших квадратов. При этом выполняется условие консервативности:
J P,dx = f |Voli|.
Vol,
Для линейной (не ТВД) схемы можно положить средние значения решения ff ) по каждой
7 " (l)
грани l ячейки равными средним значениям реконструкционного многочлена р. '.
Вычисляя поверхностные интегралы от базисных функций, получаем выражение для ^ в виде линейной комбинации значений функции /т в ячейках шаблона:
,, м
т=0
Коэффициенты
юг7т не зависят от значений решения и определяются только геометрией пространственной сетки. Для подавления паразитных осцилляций в областях быстрого изменения расчетных величин в блок реконструкции вводится так называемый ограничитель наклонов • Для расчета стационарных течений удобен ограничитель наклонов, предложенный в [27].
Окончательно средние значения по граням ячейки /1 ^, используемые в вычислениях, выражаются через средние значения (3.7) по формуле
( м
fll) = f + W, (p\l) - f ) = f + W, £ ®ilmfim - f
V m=0
(3.8)
При этом схема первого порядка соответствует выбору у, = 0, в то время как у, = 1 приводит к линейной (не ТВД) схеме второго порядка аппроксимации по пространству.
Описанная процедура реконструкции значений функции распределения на неструктурированной сетке также позволяет использовать полиномиальные представления более высокого порядка. Однако проведенные тестовые расчеты показали, что кусочно-линейное представление, соответствующее схеме второго порядка аппроксимации, обеспечивает наилучшее соотношение между точностью получаемых результатов и затратами машинного времени.
3.2.3. Нахождение потоков через грани
Вычисление численных потоков (3.5) начинается с построения для каждой ячейки Vol, вектора p, (x, y), каждая из компонент которого является реконструкционным многочленом, аппроксимирующим соответствующую компоненту вектора W в пределах ячейки. Каждая компонента p, (x, y) строится путем применения описанного выше скалярного метода реконструкции.
Затем для каждой стороны ячеек вычисляются средние значения W(l) по формуле (3.8), примененной к каждой компоненте.
В результате реконструкции для каждой грани l пространственной ячейки , получаются два
значения вектора консервативных переменных на грани. Одно значение, обозначаемое далее U-, соответствует реконструкционному многочлену ячейки Второе значение U+ получается при использовании многочлена из ячейки с номером с1 (,), прилегающей к грани l ячейки В общем
случае эти два значения не совпадают, и на грани естественным образом возникает задача о распаде разрыва [28]. Подробное описание ее аппроксимации в разностных схемах с использованием свойства инвариантности конвективного потока относительно поворота системы координат может быть найдено, например, в [29]. Пусть Т — матрица поворота системы координат в локальную систему для грани l. Тогда численный поток (3.5) имеет вид (индекс временного шага n опущен):
ф, = J Т-'F(Ти)dA = Т-f(Q-,Q+)|A|, Q± = Ти±. (3.9)
Величина Q является вектором консервативных переменных в системе координат, ориентированной по грани, функция f — так называемый решатель задачи о распаде разрыва, значения
компонент скорости вращения системы координат Vx, Vy берутся в центре грани. В выражении для вектора конвективного потока F в формуле (3.9) x-компонента скорости вращения Vx интерпретируется как проекция на нормаль к грани.
В настоящей работе используются несколько одномерных решателей задачи о распаде разрыва. Самый простой из них — поток Русанова, который определяется по формуле:
f = 2 (f (q+) + F ( Q-))-1 s (q+- Q-). (3.10)
Здесь s = max (( - Vn |, |SR - Vn |) — оценка максимальной скорости распространения волн в подвижной системе координат; SL, SR — оценки скоростей волн в задаче о распаде разрыва в неподвижной системе координат; Vn = Vxnx + Vyny + Vznz — скорость движения системы координат по нормали к грани.
Более аккуратным является метод HLL [30], модификация которого для вращающейся системы координат имеет вид:
r F-- VnU-, Sl > Vn,
F +- Vn и+, SR < Vn, (3.11)
-SL (Uhl1 -U-)-VnUhl1, SL < Vn <SR, где
U hll = SRU +- SLU + F - F + (3 12)
SR - SL
Отметим, что в (3.11), (3.12) используется выражение для конвективного потока I7 без учета вращения системы координат.
3.2.4. Решение уравнений для приращений вектора консервативных переменных
Прямое решение линейной системы (3.6) уравнений для приращений вектора консервативных переменных наталкивается на две трудности. Первая трудность состоит в вычислении матрицы Якоби SRn /SU при использовании нелинейной (ТВД) схемы второго порядка аппроксимации. Вторая трудность состоит в решении полученной системы уравнений большой размерности. В настоящей работе используется неявный метод LU-SGS без явного вычисления матриц Якоби (так называемый matrix-free метод) [31, 32]. В левой части схемы оператор дифференцирования
по пространству Rn аппроксимируется схемой первого порядка с методом решения задачи разрыва (3.10). После упрощений приращение вектора консервативных переменных U находится из следующей системы:
At
DAUi + ^КШ1AF*l -s,lAUCl(i))) = AtRn,
21 Vol, |
D = 1 + ^VTTl S s,lAl, aq = Til AU, (3.13)
2lVo1il /
AF = F (Q . (0+AQ.l (i))-F ( (i)).
Здесь sii — оценка скорости волны в задаче о распаде разрыва на грани l. После приближенной LU-SGS-факторизации получаем двухшаговую процедуру решения системы (3.13). На первом шаге обратным ходом строится промежуточное решение AU*, i = Ntot, NtOt -1, ..., 1:
At" t \
DAU*=-2\t_M 2 (Ti- 1AFii -SilAUCi(i)])ii +AtR". (3.14)
I Ml: oi (i )<i
Далее, окончательные значения вычисляются с помощью прямого хода i = 1, 2,..., Ntot:
At" / \
D AUi = Di AU* - -T-— 2 (1AF,i - Si AUCl(,))). (3.15)
2lVO4l: ol (i )>i
Число операций в численном решении системы (3.14), (3.15) линейно пропорционально числу пространственных ячеек NtOt. Реализованная неявная схема не требует существенных затрат машинной памяти и является более эффективной, чем типичная двухшаговая явная схема продвижения по времени.
Для уменьшения времени счета вычисления пространственная сетка разбивается на блоки с помощью пакета METIS [33]. С практической точки зрения наиболее предпочтительной является программная реализация метода без обмена данными между блоками на этапе решения линейной системы метода LU-SGS. При этом на границах между блоками пространственной сетки используется условие равенства нулю приращения решения. Получающийся при этом параллельный алгоритм не является эквивалентным однопроцессорному (одноблочному). Одной из задач настоящей работы является определение эффективности такого простого параллельного обобщения неявного алгоритма на многоблочные сетки.
3.2.5. Выбор шага по времени и критерии сходимости
В расчетах величина шага по времени At" вычислялась по формуле:
h
dt" = Cmin—^, (3.16)
i S"
i
где — оценка максимального (по абсолютному значению) собственного числа в ячейке; С — заданное число Куранта; И1 — характерный линейный размер ячейки. Явная схема соответствует С < 1/3. Для неявной схемы (3.4) возможно использование больших чисел Куранта С ^ 1, что ускоряет сходимость к стационарному течению.
В качестве проверки сходимости к стационарному решению во вращающейся системе координат использовалось условие малости невязки в нормах , 1^, а также установление тяги и крутящего момента. Здесь вектор невязки в ячейке 7 определяется следующим образом:
R" = 2 Ф" +S". i IVolilt i
Общая невязка для поля течения
ЕП
^ = max
R"
e" =2| R fi
Vol,.
Разностное решение считалось сошедшимся к стационарному режиму при падении обеих невязок на несколько порядков от начального значения.
3.3. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ
Описанный выше численный метод решения уравнений Эйлера во вращающейся системе координат является частью разрабатываемого пакета «Несветай-ЗД», состоящего из пяти основных частей: препроцессора, вычислительного ядра, аэродинамического решателя, кинетического решателя, а также постпроцессора.
Основной частью пакета является вычислительное ядро, разработанное В. А. Титаревым в 2007—2011 гг. [25, 26, 34, 35]. Ядро является библиотекой из четырех модулей на языке Fortran-90 и реализует базовые операции по чтению неструктурированных сеток в различных форматах, созданию информации о связности сетки, нелинейную реконструкцию функций (решения) с произвольным порядком аппроксимации, неявный метод дискретизации по времени типа LU-SGS, вывод данных для программ визуализации, а также функции обмена данными по технологии MPI. Общий размер ядра составляет примерно 10 тыс. строк кода.
Аэродинамический решатель пакета является одной из надстроек над ядром «Несветай-3Д». Другой возможной надстройкой является решатель кинетического уравнения Больцмана с модельным интегралом столкновений [26, 34]; данный решатель может применяться для численного моделирования задач механики разреженного газа.
Задача препроцессора пакета состоит в подготовке расчетной сетки, полученной из сеточного генератора, для использования на многопроцессорных ЭВМ. Постпроцессор позволяет собрать из многоблочного вывода данных одноблочный, что удобно для задачи визуализации результатов, а также для обработки данных. Отметим, что как препроцессор, так и постпроцессор используют модули ядра. Это обеспечивает автоматическое обновление всех частей пакета при обновлении ядра.
Для решения задач обтекания винта в 2012 г. была проведена модификация аэродинамического ядра для вращающейся системы координат. Полученный пакет программ позволяет рассчитывать обтекание вращающегося винта произвольной формы с применением как структурированных, так и неструктурированных сеток и является внутренним кодом ЦАГИ.
3.4. ПРИМЕРЫ РАСЧЕТОВ
Тестирование пакета проводилось на задаче обтекания вращающегося шестилопастного винта. Во всех случаях использовались нелинейная схема второго порядка аппроксимации по пространству, решатель распада разрыва HLL и неявный метод продвижения по времени. Расчеты выполнялись на суперкомпьютере «Ломоносов», установленном в НИВЦ МГУ. В качестве критерия сходимости к стационарному режиму принималось падение интегральной невязки до
величины 10-7.
3.4.1. Масштабирование кода и влияние числа Куранта
Первая серия расчетов проводилась с целью установить масштабируемость программного комплекса по числу процессоров, а также изучить влияние числа Куранта на скорость сходимости к стационарному решению. Использовалась тетраэдральная сетка из 5 млн ячеек, которая разбивалась на 32, 64, ..., 1024 блока. Так, на рис. 3 представлена поверхностная сетка, соответствующая 128 блокам. Хорошо видно сгущение сетки к кромке лопасти, а также к поворотному механизму. Полученная зависимость ускорения счета от числа ядер относительно расчета на 32 ядрах представлена на рис. 4. Видно, что программный комплекс обладает достаточно высокой эффективностью. На 1024 ядрах ускорение составляет примерно 80% от идеального. При этом число шагов по времени, требуемых для получения стационарного режима, практически не зависит от числа блоков, на которые разбита сетка.
Также проводились расчеты на гибридной сетке, состоящей из призматических элементов вблизи поверхности винта и тетраэдров в основной части расчетной области. Разбиение таких сеток с помощью METIS приводит к разбалансировке блоков по числу ячеек примерно на 6%, что несколько снижает масштабируемость метода, которая, тем не менее, остается хорошей.
I-1-1-1-1-1-1-1-1
32 280 528 776 1024
Число ядер
Рис. 4. Масштабируемость кода на сетке из 5 млн тетраэдров (расчет на 32, 128, 512 и 1024 ядрах)
Рис. 3. Поверхностная сетка для тестового расчета на 128 ядрах
О 10 000 20 000 30 000
Число шагов
Рис. 5. Влияние числа Куранта на сходимость к стационарному режиму
Кривые поведения интегральной невязки для чисел Куранта С = 0.3, 5 и 25 приведены на рис. 5 для ио = 39 м/с и числа оборотов в минуту N = 5000. Видно, что использование неявного метода продвижения и большого шага по времени приводит к существенному ускорению сходимости. Подчеркнем, что приводимые результаты относятся к полностью нелинейной схеме с включенным ограничителем наклонов. Далее в представленных расчетах принималось С = 25.
3.5. РАСЧЕТ АКУСТИЧЕСКОГО ПОЛЯ, ФВХ-ПОДХОД
Выполненные аэродинамические расчеты позволяют вычислить шум в дальнем поле, создаваемый винтом на данном режиме работы. Для получения дальнего поля используется интегральный метод, основанный на формуле Фокс Вильямса — Хоукингса (ФВХ). Эта формула представляет собой основной инструмент, с помощью которого в современных вычислительных пакетах определяется шум в дальнем поле. В своей общей формулировке данный подход позволяет учесть все источники аэродинамического шума. В качестве поверхности ФВХ, на которой задаются давление и скорость течения, в данной работе использовалась поверхность лопастей винта (рис. 6).
Рис. 6. Геометрия задачи и поверхность ФВХ /(х, Г)
Важно отметить, что использованное выражение для формулы ФВХ, полученное впервые в работе [36], в явном виде учитывает наличие равномерного среднего потока, что делает его особенно удобным при расчете шума в задачах типа «виртуальная аэродинамическая труба». Выражения для шума вытеснения и нагрузки при этом принимают следующий вид [36]: а) шум нагрузки
4пр'ь =-1
С0 г=0
Ц]П]Я^ + Ц}п ]Я, + Ц]П]Я^
■1Г
С0 3/=0
Я*(1 - М я)
дМ Я
й п-— Г
с0* /=0
дЯ
Ц]П]Я
Ц]П]Я
дт Я*(1 - Мя )3
й п + Г
б) шум вытеснения
4прТ = Г
-Г
3/=0
/=0
йп
й п +
Я*( - М я )
дМ я
й п + Г
3/=0
/=0
-дЯ
дт Я *2 (1 -М я)2 й п;
йп -
Ц]П]Я*
Я*2 (1 - МЯ )
(3.17)
йп
Я*(1 - М Я )3
дт
йп- М0 Г
0 3/=0
дт Я*2 (1 - МЯ )2
Яйп + Яй п + Яйп
й п -
+ М0 Г
/=0
дя'I Яйп
дт Я*2 (1 - Мя )2
йп-М0 Г
0->/=0
ге1
Я*(1 - М я)2
дМ я Яйп
ге1
(3.18)
- и0 Г
/=0
яйп
Я*2 (1 - Мя )2
дт Я*(1 - Мя )3 й п.
ге1
йп-
ге1
Подынтегральные значения в квадратных скобках берутся при значении времени задержки те, что в данном вычислительном комплексе было реализовано с помощью метода Advanced Time [37].
В формулах (3.17) и (3.18) введены обозначения:
R*=)j(xi -У1 )2 +Р2 (x2 -У2)2 + (хз -Уз)2 \ =-2- (-М0 +
Я = ?(о + ) Я2 = ^ = ^, д* = Х1 - у1 Х2 - у2 .^*=в2 Х3 - у3
1 я* ' 2 я* ' 3 Я* '
М я = -1 ^Д..
со
Здесь ио — скорость однородного набегающего потока; Мо — число Маха набегающего потока; п. и т — пространственная и временная координата точки поверхности ФВХ в системе отсчета источника; х. и у. — точки положения наблюдателя и источника в системе координат наблюдателя; V. — компонента скорости движения точки поверхности ФВХ в системе координат наблюдателя.
Источниковые члены Qi и Ь^ имеют вид:
Q, =Р fi (u, -v, ) + PoV, L =PfiU ( -Vj ) + (p - po^,
где и1 — скорость течения, создаваемая винтом в точке поверхности ФВХ. Вкладом квадруполь-
ной составляющей в данной задаче, по аналогии с [36], в настоящей работе пренебрегается. Суммарное звуковое поле в дальней области, таким образом, является суммой шума нагрузки (3.17) и шума вытеснения (3.18). Значения ui ир на поверхности ФВХ брались из численного расчета и не зависели от времени, так как аэродинамический расчет дает стационарное решение во вращающейся системе координат. Для данной задачи скорость точек V. самой поверхности ФВХ описывается выражениями (3.2).
Таким образом, расчет аэродинамики винта, описанный в п. 3.1 — 3.4, и расчет шума с помощью поверхности ФВХ (п. 3.5) представляют собой единый комплекс программ численного моделирования задач аэродинамики и акустики винтов, позволяющий исследовать различные конфигурации с точки зрения как шума на местности, так и шума в салоне.
4. ВАЛИДАЦИЯ И ВЕРИФИКАЦИЯ СОЗДАННОГО РАСЧЕТНОГО МЕТОДА
Для окончательного вывода о работоспособности созданного комплекса аэродинамический блок будет в этом разделе валидирован путем сравнения данных расчета с экспериментально полученными интегральными характеристиками тяги и момента. Верификация вычислительного комплекса с точки зрения акустики будет проведена путем сравнения с точным решением п. 2 для шума нагрузки и шума вытеснения, полученным в линейном приближении.
4.1. ВАЛИДАЦИЯ АЭРОДИНАМИЧЕСКОГО РАСЧЕТА ПО ТЯГЕ И КРУТЯЩЕМУ МОМЕНТУ
Для и о = 39 м/с и числа оборотов в минуту N = 5000 были проведены тестовые расчеты с размером расчетной области за винтом 2, 4 и 8 диаметров винта. Получено, что длина области 4 диаметра достаточна. Также были проведены расчеты на сетках со значительным сгущением
Рис. 7. Основная расчетная сетка из 7 млн тетраэдров
Рис. 8. Расчет аэродинамических характеристик вращающегося винта для различных параметров течения
Рис. 9. Поверхностное распределение числа Маха
вблизи винта. Общее число ячеек варьировалось от 7 до 15 млн. Получено, что разрешение сетки 7 млн достаточно для воспроизведения тяги и момента. На рис. 7 представлено сечение сетки плоскостью симметрии, проходящей через ось винта. Вычисления проводились на 128 ядрах.
На рис. 8 представлены результаты расчета тяги и крутящего момента для фиксированной скорости вращения винта 5000 об./мин и различных значений скорости набегающего потока в сравнении с экспериментальными данными. Следует отметить, что расчеты проводились на основе уравнений Эйлера без учетов эффектов вязкости, теплопроводности и турбулентности. Отклонения вычисленных значений тяги и крутящего момента от экспериментальных данных максимальны при малых скоростях набегающего потока и составляют около 5% от пиковых значений данных величин. Для больших скоростей потока различие составляет менее 1%. С учетом указанных выше факторов такое согласие с экспериментом можно считать вполне приемлемым.
В заключение раздела рис. 9 иллюстрирует распределение числа Маха по поверхности винта для скорости набегающего потока ио = 39 м/с.
4.2. ВЕРИФИКАЦИЯ РАСЧЕТНОГО МЕТОДА ПО ЛИНЕЙНОМУ ПРИБЛИЖЕНИЮ
Для верификации программы расчета шума было проведено сравнение расчета шума винта в дальнем поле на основе подхода ФВХ (формулы (3.17), (3.18)) с результатами п. 2. Для расчета шума винта по формуле (2.31) в 50 сечениях лопасти вдоль радиуса (рис. 10) брались геометрические параметры (хорда, максимальная толщина лопасти в данном сечении и функция распределения толщины лопасти вдоль хорды) и рассчитанные нагрузки (функция распределения нагрузок вдоль хорды, суммарное значение нагрузок в данном сечении).
Рис. 10. Вид лопасти винта в плане (черные линии соответствуют сечениям лопасти)
90 85 80 75 70 65 60 55 50
55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 Поток i-¡> —ФВХ -«-Линейное приближение
Рис. 11. Сравнение направленности амплитуды тональной гармоники на основной частоте (500 Гц) для численного расчета (•) и линейного решения (▲) (угол 90° соответствует плоскости вращения винта)
Для сравнения подходов использовались амплитуды тональных гармоник на основной частоте винта (500 Гц). Результаты акустических расчетов приведены на рис. 11. Легко видеть, что результаты расчетов по обеим методикам показывают хорошее согласование результатов. Расхождение между линейным решением и ФВХ составляет меньше 1 дБ в плоскости вращения винта и увеличивается вверх по потоку в связи с ухудшением точности верифицирующего решения (2.31) при малых углах.
5. ЗАКЛЮЧЕНИЕ
Разработанный вычислительный комплекс программ позволяет создать среду моделирования, включающую в себя все этапы решения прикладных задач: получение геометрической модели объекта, построение расчетной сетки в физическом пространстве, расчет обтекания с использованием суперкомпьютеров, обработку результатов в акустическом постпроцессоре и выдачу рекомендаций разработчикам. Получено хорошее согласование с экспериментом по аэродинамическим характеристикам для вращающегося винта реальной формы. Валидация расчетов по имеющимся экспериментальным данным тяги и момента показала высокую эффективность метода. Для верификации акустических расчетов дальнего поля проведено сравнение с точным решением линейного приближения, в котором краевая задача для давления сведена к неоднородному конвективному волновому уравнению с известным поверхностным распределением монопольных и дипольных источников. Сравнение численного метода с точным решением показало хорошую точность разрабатываемого метода и с точки зрения акустических расчетов.
Авторы выражают благодарность А. Ф. Чевагину за поддержку этих работ в ЦАГИ, а также А. В. Лысенкову и А. А. Савельеву (ЦАГИ) за помощь в обработке геометрии винта и полезные обсуждения результатов аэродинамических расчетов. Авторы выражают признательность рецензенту за высказанные замечания, работа над которыми существенно улучшила статью.
ПРИЛОЖЕНИЕ А. ЛОПАСТЬ С ОТГИБОМ
Пусть форма поверхности невозмущенной лопасти будет Е = Е0 (р), где Е0 (р) — произвольная функция. Уравнение верхней и нижней поверхности имеет вид
Е+=Ео + /+( р), Е-=Ео +¥~(У, Р).
Ковариантные компоненты внешней нормали к верхней и нижней поверхности с использованием метрического тензора (2.4) будут
n±=±
0i
0
W
df ±/dí df ±/dp
(А.1)
где
% =
Г1 ^ о
Ч^о )
> "о
п0 = яипо =
о ]
Ы = ё']п0г п0 ] = 1
Тогда граничное условие (2.9) примет вид:
(поуф)
дф дп
и
я11 _ ео я
я12 _ ео я
я13 _ео
и2 °1
и2 и Л Р)
дк
Л л дУ
13 ^ 23
(А.2)
(А.3)
(А.4)
е' _ ио ео ^о тт2 и2 Р
Для перехода от скачка разований получим
отсюда
дф дп
др дп
используем формулу (2.13). После некоторых преоб-
(У ) = Ро (идУ)( поУф),
др дпо
Р и2 д2 к
РоиЛ
(
ио ео ео т т2 и2 Р
(А.5)
Отличие формулы (А.5) от аналогичной формулы (2.15) для лопасти без отгиба связано с необходимостью нормировки по в соответствии с (А.3). Плотность источников, соответствующую (А.5), определим с помощью (2.17):
(прУр) 5
4 =11 II 5.
Можно показать, что 5^ = ||по || 5 (2, _ ео (г)). Отсюда
д2г
д = (поУр)5 ( _ ео (г)) = -РоиЛ — 5 ( _ ео (Р)),
дУ
(А.6)
т. е., несмотря на различие в формулах для скачка давления (А.5) и (2.15), источники в волновом уравнении тем не менее оказываются одинаковыми.
Для плотности источников шума нагрузки из формулы (2.17) получается:
(поУ)([рК) = п
■у, [ р (у, г)ы 5(-ео (р))] .
(А.7)
/'о|| ||по|' .....
Формула (А.7) после раскрытия имеет очень громоздкий вид, так как для нормали по приходится использовать контравариантные координаты п'о (А.2). Необходимо обратить внимание,
1
1
п
о
что формула (А.7) для источников в случае шума нагрузки состоит не только из членов, пропорциональных 5'( —^0 ), как в случае плоской задачи, или лопасти без отгиба (2.27). Она включает
также члены, пропорциональные 5 ( — ^0 ).
ПРИЛОЖЕНИЕ Б. СРАВНЕНИЕ ПОЛУЧЕННОГО РЕШЕНИЯ С ФОРМУЛОЙ ХАНСОНА
Сравним полученные результаты для дальнего поля винта с известными результатами работы Хансона [12]. В этой работе использовалась система координат, связанная с Землей. В связи с этим приведем результаты настоящей работы к неподвижной системе координат.
Пусть источник, находящийся в точке А (рис. 12), создает звуковой сигнал. Через некоторый интервал времени t0) источник, пройдя расстояние , окажется в точке B, а сигнал за это время достигнет точки наблюдения C, находящейся на расстоянии от точки A. Координаты r, 9 использовались выше при рассмотрении звукового поля в движущейся системе. Они представляют собой сферические координаты с центром в точке B. В системе координат, связанной с Землей, как правило, используются координаты Г, 9, представляющие собой сферические координаты с центром в точке A. Отметим, что координаты r, 9 во всем пространстве отсчиты-
Рис. 12. Системы координат
ваются от общего центра в точке В, где в рассматриваемый момент времени находится источник. В то же время координаты г, 0 отсчитываются от различных центров для разных точек наблюдения, поскольку к этим точкам одновременно приходит сигнал, созданный движущимся источником на различных участках всего пути. Для неподвижной точки наблюдения обе системы г, 0
и Г, 0 являются, очевидно, функциями времени. Координаты г, 0 иг, 0 связаны между собой соотношениями:
1
г = a • r\ 1 —
1
— M2 sin2 9,
a=
1 — M0 cos 9 1 — M2
M0 cos 9 — M2 sin2 9
(Б1)
Рассмотрим дальнее звуковое поле винта в неподвижной системе координат. Переходя к переменным г, 0, получим, что п-я гармоника шума винта в области дальнего поля
pn (r, t) = —
jng—iat + ikr
4nr (1 — M0 cos 0)
^ (9)
fmp
(Б.2)
где величина
^n (9)
для шума вытеснения
^ (0) = \ак (аМлкп ) (Ькпр')р'
(Б.3)
и для шума нагрузки
^ (0 ) = _{ ,пи Л
и оР'
1 _ а
и Л
ро (аМлкп) (кУЖ
(Б .4)
1
а = -
1 _ Мо 008 0
Ь = •
8Ш 0
1 _ Мо 008 0
Как уже отмечалось, координаты г, 0 для неподвижной точки наблюдения являются функциями времени. Однако для достаточно малого промежутка времени наблюдения г (такого, за который источник проходит путь, много меньший, чем расстояние от него до точки наблюдения) Г (г ) и 0 (г) можно разложить в ряд
:(г) = Го _
008 0о
1 _ Мо 008 0о
-уг-
( Б .5)
где через го и 0о обозначаются г (г) и 0(г) при г = о.
Именно координаты Го и 0о используются в [12] при вычислении звукового поля винта в системе координат, связанной с Землей. Подставляя (Б.5) в (Б.2) и оставляя только члены, неубывающие при Го получим:
ехр
рп ( г) = _-
гапг
1 _ М П 008 0о
,кпГо
4пго (1 _ Мо 008 0о)
Р (0о )
,пе,пф.
(Б.6)
В этом выражении на первый взгляд кажется, что фазовая скорость не совпадает со скоростью звука. Однако это не так. Фазовая скорость не может быть рассчитана здесь исходя из отношения коэффициентов при временном и пространственном слагаемых под экспонентой, поскольку координаты г, 0 отсчитываются от различных центров для разных точек наблюдения, и к этим точкам одновременно приходит сигнал, созданный движущимся источником на различных участках всего пути.
Аналогичный результат, полученный в работе [12], несколько отличается от (Б.6):
(
ехр
рп( г) = _"
/<впг
,кпГо
\
1 _ Мо 008 0о 1 _ Мо 008 0
о )
4пго (1 _ Мо 008 0о )
Р(0о )
,Пе,Пф.
( Б .7)
Действительно, сравнение (Б.7) и (Б.6) показывает, что вместо набега фазы от источника до
точки наблюдения, равного кпго, в [12] получен набег фазы
кпГо
1 _Мо 008 0о
что при немалых М
может приводить к сильным фазовым искажениям в звуковом поле. Дело в том, что разложение плотности источников я по малому параметру г '/г , использованное в [12] для случая дальнего поля, применимо только в системе координат, движущейся вместе с источником, так как в этой системе источник занимает ограниченную область пространства. Использование такого разложения в [12] в неподвижной системе координат, где область, занимаемая источником, за все время движения неограниченна, и привело к указанному выше выражению.
ЛИТЕРАТУРА
1. Гутин Л. Я. О звуковом поле вращающегося винта // ЖТФ. 1936. 6(5), с. 899—909.
2. Непомнящий Е. А. Исследование и расчет звука воздушного винта // Труды ЦИАМ, № 39. — М.: Оборонгиз, 1941, с. 78.
3. Блохинцев Д. И. Акустика неоднородной движущейся среды. — М.: Наука,
1981.
4. Mikkelson D. C., Mitchell F. G., Bober L. J. Summary of recent NASA propeller research // NASA TM 84-32344. 1984, p. 1—42.
5. Ганабов В. И., Мунин А. Г. О расчете шума вращения одиночного винта с лопастями произвольной формы // Ученые записки ЦАГИ. 1989. Т. XX, № 5, с. 43—52.
6. Мунин А. Г. Авиационная акустика. — М.: Машиностроение, 1986.
7. Зленко Н. А., Кедров А. В., Кишалов А. Н. Оптимальное аэроакустическое проектирование воздушного винта // Ученые записки ЦАГИ. 2011. Т. XLII, № 6, с. 92—103.
8. Tam C. K. W., Salikuddin M. Weakly nonlinear acoustic and shock-wave theory of noise advanced high-speed turbopropeller // J. Fluid Mech. 1986. V. 164, p. 127—154.
9. Farassat F. Linear acoustic formulas for calculation of rotating blade noise // AIAA J. 1981. V. 19, N 9, p. 1122—1130.
10. Magliozzi B., Hanson D. B., Amiet R. K. Aeroacoustics of flight vehicles. Theory and Practice / Ed. by H. H. Hubbard // Acoust. Society of Amer., 1991. V. 1, p. 1 —64.
11. Hanson D. B. The aeroacoustics of advanced turbopropellers // Mechanics of sound generation in flows. — Springer, 1979, p. 282—293.
12. H a n s o n D. B. Helicoidal surface theory for harmonic noise of propellers in the far field // AIAA J. 1980. V. 18, N 10, p. 1213 — 1264.
13. Hanson D. B. Compressible helicoidal surface theory for propeller aerodynamics and noise // AIAA J. 1983. V. 21, p. 881—889.
14. Ffowcs Williams J. E., Hawkings D. L. Sound generation by turbulence and surfaces in arbitrary motion // Phil. Trans. Royal Soc., Series A. 1969. V. 264, p. 321 —342.
15. Bosnyakov S., Kursakov I., Lysenkov A., Matyash S., Mikhai-lov S.,Vlasenko V.,Quest J. Computational tools for supporting the testing of civil aircraft configurations in wind tunnels // Progress in Aerospace Sciences. 2008. V. 44, p. 67—120.
16. Yin J., Stuermer A., Aversano M. Coupled URANS and FW-H analysis of installed pusher propeller aircraft configurations // AIAA Paper 2009-3332.
17. Deconinck T., Capron A., Barbieux V., Hirsch C., Ghorbaniasl G. Sensitivity study on computational parameters for the prediction of noise generated by counter-rotating open rotors // AIAA Paper 2011-2765.
18. Barbarino M., Casalino D. Hybrid analytical/numerical prediction of propeller broadband noise in the time domain // Int. J. Aeroacoustics. 2012. V. 11, N 2, p. 157—176.
19. Kopiev V. F., Maslov A. A., Chernyshev S. A. Prop-fan sound field shielding by the fuselage boundary кует // DGLR/AIAA Paper 92-02-068. 1992, p. 5.
20. Владимиров В. С. Уравнения математической физики. — М.: Наука, 1976,
528 с.
21. H i r s h C. Numerical computation of internal and external flows // Elsevier, 2007.
22. Biedronand R. T., Thomas J. L. Recent enhancements to the FUN3D flow solver for moving-mesh applications // AIAA Paper 2009-1360.
23. Dumbser M., Käser M. Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems // J. Comput. Phys. 2007. V. 221, N 2, p. 693—723.
24. Dumbser M., Käser M., Titarev V. A., Toro E. F. Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems // J. Comput. Phys. 2007. V. 226, p. 204—243.
25. Tsoutsanis P., Titarev V. A., Drikakis D. WENO schemes on arbitrary mixed-element unstructured meshes in three space dimensions // J. Comput. Phys. 2010. V. 230, p. 1585—1601.
26. Titarev V. A. Efficient deterministic modeling of three-dimensional rarefied gas flows // Commun. Comput. Phys. 2012. V. 12, N 1, p. 161 — 192.
27. V e n k a t a k r i s h n a n V. On the accuracy of limiters and convergence to steady-state solutions // AIAA Paper 93-0880, 1993.
28. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики // Мат. сборник. 1959. Т. 47, № 89, с. 271 —306.
29. Toro E. F. Riemann solvers and numerical methods for fluid dynamics. — SpringerVerlag, 2009.
30. Harten A., Lax P. D., van Leer B. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws // SIAM Review. 1983. V. 25, N 1, p. 35—61.
31. Men'shov I. S., Nakamura Y. An implicit advection upwind splitting scheme for hypersonic air flows in thermochemical nonequilibrium // A collection of technical papers of 6th Int. Symp. on CFD. 1995. V. 2, p. 815.
32. Men'shov I. S., Nakamura Y. On implicit Godunov's method with exactly linearized numerical flux // Computers and Fluids. 2000. V. 29, N 6, p. 595—616.
33. Kary p i s G., Kumar V. Multilevel k-way partitioning scheme for irregular graphs // J. Parallel Distr. Com. 1998. V. 48, p. 96—129.
34. Титарев В. А. Неявный численный метод расчета пространственных течений разреженного газа на неструктурированных сетках // ЖВМиМФ. 2010. Т. 50, № 10, с. 1811—1826.
35. Titarev V. A., Drikakis D. Uniformly high-order schemes on arbitrary unstructured meshes for advection-diffusion equations // Computers and Fluids. 2011. V. 46, N 1, p. 467—471. 10th ICFD Conference Series on Numerical Methods for Fluid Dynamics (ICFD 2010).
36. Najafi-Yazdi A., Bres G., Mongeau L. An acoustic analogy formulation for moving sources in uniformly moving media // Proc. Royal Soc., A. 2011. V. 467, N 2125, p. 144—165.
37. C a s a l i n o D. An advanced time approach for acoustic analogy predictions // J. Sound Vibration. 2003. V. 261, p. 583—612.
Рукопись поступила 15/II2012 г.
