Анализ звукового удара с использованием полей возмущений, рассчитанных по нелинейной теории Текст научной статьи по специальности «Физика»

Том X Ь

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2 0 09

№ 4

УДК 534.83:629.7.015.016.54

АНАЛИЗ ЗВУКОВОГО УДАРА С ИСПОЛЬЗОВАНИЕМ ПОЛЕЙ ВОЗМУЩЕНИЙ, РАССЧИТАННЫХ ПО НЕЛИНЕЙНОЙ ТЕОРИИ

ЧО КЮЧУЛ

Разработана методика построения эпюры звукового удара на местности, основанная на расчете возмущенных полей течения в ближнем поле самолета при численном решении полных уравнений Навье — Стокса. В основе алгоритма лежит методология расчета звукового удара, развитая в ЦАГИ Жилиным Ю. Л. и его учениками. Применение нелинейных методов расчета возмущенных полей течения позволяет учесть в той или иной степени такие факторы, влияющие на интенсивность звукового удара, как затупление носовой части, вязкость, струи двигателей и т. д. В качестве примера показано влияние струи силовой установки самолета Ту-144 на эпюру звукового удара, дано сравнение с экспериментом.

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

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

Рис. 1. Процесс вычисления интенсивности звукового удара на земле

фюзеляж — крыло [3]. В основе метода лежит вычисление некоторой функции Е, названной именем Уизема. Исходными данными метода является распределение площади поперечных сечений и подъемной силы по длине компоновки самолета. В правой части рис. 1 показан ход вычисления звукового удара. Возмущенное поле течения около реального объекта заменяется полем от некоторого эквивалентного тела вращения. Для возмущений, распространяющихся непосредственно под самолетом, площадь поперечных сечений эквивалентного тела вращения состоит из двух слагаемых — геометрического 8а и аэродинамического :

Распределение эквивалентной площади 8Л (х) без подъемной силы определяется по сечениям самолета, наклоненным по характеристическим линиям Маха. Полученные сечения проектируются на плоскость х = 0, где вычисляются их площади. Распределение площади эквивалентного тела в зависимости от распределения подъемной силы по длине самолета равно:

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

к — отношение теплоемкостей. Функция Е, пропорциональная избыточному давлению, определяется как производная потенциала от распределения источников вдоль оси эквивалентного тела вращения Ф:

^ (х) = 8л (х) + ^т (х).

(1)

(2)

где в = -у/М^ -1 ; Мте — число Маха набегающего потока; q — скоростной напор; Т(х) — подъем-

кМ

(3)

Е (п) = з^ф(п)

ОП

(4)

где п = п/^ — расстояние вдоль оси тела, отнесенное к характерной длине. Функция Е (п) получается дифференцированием функции Ф(п) по п. На больших расстояниях от самолета Ф(п) вычисляется при помощи асимптотического преобразования [1, 2]:

(5)

Зависимость избыточного давления Ар от времени ^ определяется соотношениями [4]:

(6)

Р_2 - Ъ = 1, П2 - П1 к

_ _ Л2 _

( + Ъ)((2 -П ) = 21 ^, (9)

Л1

где индексы 1 и 2 при функциях Ъ означают, что соответствующая величина вычисляется на характеристической поверхности П = П1 или П=П2- Коэффициенты к, к и ^2 [4], определяющие затухание возмущения в слоистой атмосфере, зависят только от состояния атмосферы и режима полета самолета и не зависят от формы самолета. Уравнения (6), (7) описывают затухание возмущений на пути от летального аппарата до земли, а (8), (9) — процесс образования и перемещения ударных волн в пакете возмущений. Уравнение (9) является следствием того, что скорость распространения ударной волны относительно частиц воздуха равна среднему значению скорости звука непосредственно до и после ударной волны. Учет зависимости скорости распространения волны от ее интенсивности приводит к более быстрому затуханию возмущений, чем в геометрической акустике. С физической точки зрения, более быстрое затухание возмущений объясняется диссипацией энергии в ударных волнах.

Система уравнений (6) — (9) может быть решена численно. Алгоритм решения системы изложен в работе [5]. Данную методику можно применять не только для определения эпюры звукового удара на местности, но и для получения избыточного давления волны на любом расстоянии от самолета. Это обстоятельство позволяет получить исходные данные для программы 2БРИУЯи8 [6], которая рассчитывает распространение звукового удара с учетом фокусировки звукового луча, отражения волны от верхних слоев атмосферы (вторичный звуковой удар), влияния атмосферной абсорбции и дисперсии.

Метод вычисления распределения площади эквивалентного тела вращения по интегральным соотношениям. Для определения эквивалентного тела вращения можно воспользоваться интегральным соотношением, полученным в работе [7] для слабо возмущенного сверхзвукового потока. Этот метод целесообразно применять для построения эквивалентного тела вращения компоновок самолетов с внутренними протоками и струями двигателя, а также для учета толщины вытеснения пограничного слоя. На рис. 2 показан контрольный объем сверхзвукового обтекания компоновки крыло — фюзеляж. Поверхность является границей возмущенного и невозмущен-

ного потоков (конус Маха), плоскость «2, параллельная плоскости хо^, находится под телом, но не касается его. Плоскости «3 и «4 определяются огибающими обратных конусов Маха, исходящих из отрезка АВ. Если плоскость «2 и соответственно отрезок АВ перпендикулярны оси у, распределение площади « соответствует эквивалентному телу вращения, возмущения от которого

Рис. 2. Контрольный объем при сверхзвуковом обтекании компоновки крыло — фюзеляж

распространяются непосредственно под самолетом. В этом случае рассчитать распределение производной площади можно по формуле, полученной в [7]:

dSea (х) 1 в 2В Бг 2 В

1 =-7Т | = -7Т 1 “* = ~7Т 1и*- (10)

А а ^ а

где и, V — компоненты возмущенной скорости вдоль осей х, у.

В рамках линейной теории плоскость £{, на которой вычисляется интеграл (10), может находиться на любом расстоянии от тела, в том числе и на очень малом. На рис. 3. показана эффективная длина самолета и положение плоскости под ним при полете под углом атаки а. В нелинейной сверхзвуковой теории угол скачка уплотнения ц меняется вдоль длины самолета, и в конце самолета угол принимает значение рі, которое не равно ц, что необходимо учитывать при определении эффективной длины а - Ь.

У

и

Эффективная длина Рис. 3. Эффективная длина и зона плоскости ^

Метод вычисления интенсивности звукового удара в данной работе. На рис. 4. показана схема вычисления интенсивности звукового удара. Сначала строится математическая модель поверхности, повторяющей реальную компоновку с необходимой точностью. Затем компоновка заключается в расчетную область, где строятся поверхностная и объемная неструктурированные сетки. Для точного вычисления возмущения в плоскости , где вычисляются начальные данные для задачи звукового удара, необходимо сгущать сетку не только на поверхности компоновки, но и в зоне между компоновкой и плоскостью . Также нужно создать специальную структурированную сетку около поверхности компоновки для правильного учета вязкости и построения профиля скоростей в пограничном слое. На следующем этапе формируются граничные условия на внешней поверхности расчетной области, на поверхности компоновки, на входе в воздухозаборник и на срезе сопла двигателя.

Рис. 4. Схема вычисления интенсивности звукового удара (ЗУ) в данной работе

Для интегрирования системы уравнений Навье — Стокса используется метод конечных объемов [8]. Расчетная область разделяется сеткой на контрольные объемы. Интегральные уравнения для каждого контрольного объема представляются системой алгебраических уравнений, которые решаются итерационно [8]. Турбулентность моделируется статистическими турбулентными моделями по осредненному числу Рейнольдса [9].

Интеграл (10) от возмущенной осевой скорости в плоскости 5*2 дает распределение производной площади поперечных сечений эквивалентного тела вращения, по которой определяется функция Уизема. Затем определяются коэффициенты к, к^ ^ и решается система уравнений

(6) — (9), что дает либо эпюру звукового удара на земле, либо исходные данные для программы гЕРНУЯШ.

Определение ближнего поля с помощью эквивалентного тела вращения. При расчете возмущенного поля давления в сверхзвуковом потоке около компоновки было обнаружено, что результат сильно зависит от параметров сетки между компоновкой и расчетной плоскостью. На рис. 5, а показан результат расчета поля давления при М^ = 2 и угле атаки, равном нулю, в сверхзвуковом потоке под телом вращения, состоящим из конуса и цилиндра, общей длиной Ь = 10.16 см, с радиусом миделя Я = 0.48 см.

Расчетная плоскость находится под телом на расстоянии у = -18.9 см. Сравнение численных и экспериментальных результатов [10] показано на рис. 5, б. Непосредственное вычисление возмущенного относительного давления показало значительное отличие от эксперимента. На рисунке сплошная линия обозначает результат, который был получен по разработанной методике, с использованием возмущенного давления в ближнем поле под телом (у = -0.51 см). Этот пример показывает, что комбинированный линейно-нелинейный подход к определению возмущенного поля на значительном расстоянии от тела позволяет использовать конечно-разностные методы со сгущением сетки только вблизи тела (у = -0.51 см).

о Эксперимент [10] ТТпетмпй ПЯГДІЯТ

о* Л_ Блше «нее поле и теори. я звукового удара

\

/с )

\ Г\ \\ к -О..

^ и - й\ о

о\ ог

\ N

6) -5 0 5 10 15 20 1, см

Рис. 5. Распространение возмущений в плоскости симметрии (а) и изменение давления под телом вращения (б) (у = -18.19 см)

Расчет эпюры звукового удара самолета Ту-144. Компоновка самолета Ту-144 без силовой установки. Составляющие распределения площади поперечных сечений эквивалентного тела вращения, соответствующего крейсерскому полету компоновки самолета Ту-144, состоящей из крыла, фюзеляжа и вертикального оперения показаны на рис. 6. Для этого расчета были использованы следующие основные параметры: число Маха — 1.98, длина самолета — 64.45 м,

масса — 157 т, высота полета — 16100 м, коэффициент подъемной силы (еу) — 0.1095, угол атаки — 5.35°.

Эквивалентное тело было определено по распределению площадей поперечных сечений и подъемной силы вдоль оси симметрии компоновки. На рис. 7, а показано возмущенное поле давления под самолетом Ту-144 в плоскости ^ Максимальное избыточное давление (зона Р1) в плоскости $2 индуцирует сверхзвуковая часть передней кромки консоли крыла самолета. Минимальное избыточное давление (зона Р2) индуцируется волнами, исходящими от задней концевой кромки крыла.

Максимальное и минимальное избыточное давление определяет максимальный и минимальный пик в распределении производной площади эквивалентного тела вращения $^.

На рис. 7, б показано сравнение эпюр производной площади эквивалентного тела вращения, определенных непосредственно по распределению площадей и подъемной силы в сечениях самолета плоскостями Маха (классический метод) и по возмущенному полю в плоскости $2 . Расхождение в производной площади эквивалентного тела вращения, построенной двумя способами, в зоне х = 0 ~ 0.25 объясняется нелинейным поведением ударных волн между носком фюзеляжа и плоскостью интегрирования $2, а в зоне х = 0.65 — 1 (от задней кромки крыла до конца компоновки) аналогично влияет нелинейное поведение волн и концевых вихрей. При полете под углом атаки а (рис. 8) расстояние между строительной плоскостью самолета и плоскостью $2 уменьшается. В зоне носка фюзеляжа ударные волны успевают соединиться до того, как они достигнут плоскости $2. Этот эффект объясняет отличие эпюры в зоне х = 0 — 0.25. При интегрировании возмущения между сечениями 1 и 2 в плоскости $2 возникает дополнительная зона, в которой содержатся возмущения, индуцированные ударными волнами и вихрями, сбегающими с конца крыла. Это объясняет отличие эпюр в зоне х = 0.65 — 1.

самолета Ту-144

Рис. 7. Поле возмущенного давления в плоскости $2 под самолетом Ту-144 (а)

и сравнение производной площади поперечного сечения эквивалентного тела

вращения (б)

Рис. 8. Распространение возмущения от компоновки к плоскости $2 при полете

с углом атаки а:

И — расстояние между носовой частью самолета и плоскостью $2 по У, 1, 2 — плоскости Маха на задней кромке крыла и за крылом

Интегрирование возмущенной осевой скорости на плоскости $2 позволяет определить функцию Е. Для оценки распространения возмущения в атмосфере находятся коэффициенты к, кі, ^2. На рис. 9 расчетная интенсивность звукового удара сравнивается с экспериментом [12], полученным для параметров полета, указанных ранее. В табл. 1 приведены величины максимального и минимального перепада давления, полученные по различным методикам. Результат, определенный по классическому методу (метод 2), ближе к результату метода 3, в котором обтекание компоновки и возмущенное поле в плоскости $2 получено по линейной теории [13]. Максимальные избыточные давления Артах метода 2 отличаются от эксперимента на 10.2 Па, а метода 3 — на 11.7 Па. Минимальные избыточные давления Аршіп отличаются на 21.3 Па и 25.3 Па для методов 2 и 3 соответственно.

№ Метод ^Ршах, Па ^Рт!^ Па

1 Эксперимент [12] 107 -80.7

2 Классический метод (по формуле (1)) + к, кх, к2 117.2 -102

3 Линейная теория (возмущение) [13] + 32 + к, кх, к2 118.7 -106

4 Нелинейная теория (возмущение) + 32 + к, кх, к2 113.3 -99.6

5 Нелинейная теория (возмущение) + 32 + 7ЕРНУЯи8 101.9 -95.8

Максимальное избыточное давление Артах, полученное по возмущениям, вычисленным по нелинейной теории (метод 4), меньше результата метода 3 на 5.4 Па, а Арт;п меньше на 6.4 Па; Артах по методу 4 больше, чем по методу 5, на 11.4 Па, а Арт;п — на 3.8 Па. Максимальное избыточное давление Артах по программе 2ЕРНУЯи8 меньше, чем по эксперименту, на 7 Па, а минимальное избыточное давление Арт;п меньше на 15.1 Па. Меньшие значения Артах и Арт;п по программе 2ЕРНУЯи8 можно объяснить влиянием абсорбции и дисперсии звуковых волн в атмосфере. Результаты нелинейной теории для Артах и Арт1п (методы 4, 5) в целом ближе к эксперименту, чем результаты линейной теории (методы 2, 3), хотя отличие в минимальном избыточном давлении остается значительным.

Ар, Па

I % 1 1 3 4 5 О Эксперимент [ 12] Классический метод + к, к^ к2 Линейная теория (возмущение) [13] + 52 + к,кгкг_ Нелинейная теория (возмущение) + 52 + к, кГ к2 — • — Нелинейная теория (возмущение) + 52 + 2ЕРНУЁ.и8

1 ОІ&Ч

1 V

°о < о . о 3 г1- о 6

ч.

N ’ІЧ

-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Г, с

Рис. 9. Сравнение эпюры избыточного давления в А-образной волне

Компоновка самолета Ту-144 с силовой установкой. При моделировании силовой установки можно изучить влияние струй двигателей на поле возмущения «2 и эпюру звукового удара. В данной работе для моделирования силовой установки использовались следующие параметры: площадь сопла «с — 1.35 м ; радиус сопла Я — 0.656 м; температура на срезе сопла — 670 К; скорость выходящей струи — 1160 м/с; расход воздуха — 83 кг/с; давление на срезе сопла равно давлению набегающего потока; число Маха полета М„ = 1.98; высота полета — 16.1 км [11]. Кроме этих параметров необходимо было задать площадь входа воздухозаборника, длину силовой установки, угол установки двигателей, положение силовой установки относительно крыла и т. д. На рис. 10 показана форма, положение и размеры силовой установки, использованные в данной работе.

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

Рис. 10. Схема компоновки самолета Ту-144 с силовой установкой [11]

Рис. 11. Влияние силовой установки на производную площади эквивалентного тела вращения

и функцию Уизема:

а — поле давления в плоскости ^ б — линии тока около фюзеляжа, крыла и сопла; в — распределение производной площади эквивалентного тела вращения 5^; г — распределение функции ¥

струй и их форма меняются. Сравнение производной площади эквивалентного тела вращения без силовой установки и с силовой установкой показано на рис. 11, в. Отличие площади эквивалентного тела вращения в зоне х = 0~0.3 определяется углом атаки и расстоянием к между носком фюзеляжа и плоскостью 52. В табл. 2 показаны параметры для определения положения плоскости 52. При практически одинаковой величине коэффициента подъемной силы углы атаки компоновок отличаются. Расстояние к от плоскости 52 до компоновки с силовой установкой больше, чем для компоновки без силовой установки, так как следует избегать пересечения струями двигателей плоскости 52.

Компоновка Угол атаки, град к, м

Без силовой установки 5.35 6.3 0.1095

С силовой установкой 5.23 9.3 0.1096

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

струй эпюра производной площади 5^ значительно отличается от эпюры компоновки без силовой установки. В окрестности точки х = 0.65 заметно влияние зоны давления р3, которое зависит от параметров струй, истекающих из сопла. На рис. 11, г показано влияние двигательной установки на функцию Е. Сравнение соответствующих эпюр звукового удара на земле показано на рис. 12.

Ар, Па 120

В» о Эксперимент [12] Компоновка без силовой установки

X Компоно вка с сил овои уста -ювкои

о

Ч

ц

80 40 0

-40 -80 -120

-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Г, с

Рис. 12. Влияние двигательной установки на эпюру звукового удара на земле

В результате для компоновки самолета с силовой установкой максимальное избыточное давление А^шах = 112.4 Па, а минимальное избыточное давление ДрШщ = -95.4 Па. Эти результаты отличаются от компоновки без двигателей на 0.9 Па для Арщ^ и на 4.2 Па — для АрШщ . Учет влияния струй сопла снижает величину минимального избыточного давления эпюры звукового удара, однако отличие расчетного значения от экспериментального остается. Аpmлx компоновки с силовой установкой несколько меньше, чем у компоновки без силовой установки, хотя в распределении функции Е максимум для компоновки с двигателями больше, чем для компоновки без них. Для объяснения этого целесообразно рассмотреть формулы (8), (9). Геометрическая интерпретация процесса затухания амплитуды волны в атмосфере иллюстрируется на рис. 13,

где показано распределение функции Р. Интенсивность затухания звуковой волны определяется наклоном линии 1^, соединяющей точки кривой функции Е таким образом, что площади ^, 5^2 совпадают. При учете влияния двигателей функция Е приобретает как отрицательные, так и положительные приращения (пунктирная линия), но наклон линии ^ практически не меняется. Поэтому в данном примере двигатели не оказывают существен-

Рис. 13. Геометрическая интерпретация процесса н°ГО влияния на затухание амплигуды звук°-

формирования звуковой волны вой волны компоновки.

Выводы. В работе изложена методика расчета эпюры звукового удара на местности компоновки сверхзвукового самолета при наличии двигателей. Приведены результаты расчета без силовой установки и с силовой установкой с использованием конечно-разностного метода расчета уравнений Навье — Стокса. Избыточные давления без силовой установки и с силовой установкой, полученные по нелинейной теории, ближе к экспериментальным значениям, чем результаты, полученные по возмущениям, вычисленным по линейной теории. Под влиянием струй двигателей расчетное избыточное давление в эпюре звуковой волны уменьшается по сравнению с компоновкой без двигателей.

Автор выражает глубокую признательность В. В. Коваленко за ценные замечания и помощь в оформлении статьи.

ЛИТЕРАТУРА

1. Whitham G. B. The behavior of a supersonic flow past a body of revolution far from the axis // Proc. Roy. Soc. 1950. V. 201, N 1064.

2. Whitham G. B. The flow pattern of a supersonic projectile // Communications on Pure & Applied Math. 1952. V. 5, N 3.

3. Walkden F. The shock pattern of a wing-body combination far from the flight path //

Aeronautical Qtr. 1958. V. 9. Pt 2.

4. ЖилинЮ. Л. О звуковом ударе // Ученые записки ЦАГИ. 1971. Т. II, № 3.

5. ЖилинЮ. Л., ЧернышевС. Л. Алгоритм построения эпюры избыточного давления при звуковом ударе // Труды ЦАГИ. 1981, вып. 2110.

6. Robinson L. D. Sonic boom propagation through an inhomogeneous, windy atmosphere: Ph.D. dissertation. — University of Texas, 1991.

7. ЖилинЮ. Л., КоваленкоВ. В. О связывании ближнего и дальнего полей в задаче о звуковом ударе // Ученые записки ЦАГИ. 1998. Т. XXIX, № 3 — 4.

8. ANSYS, Inc., ANSYS CFX-solver, Release 10.0, 2005: Theory, Discretisation and Solution Theory.

9. ANSYS, Inc., ANSYS CFX-solver, Release 10.0, 2005: Theory, Turbulence and Wall Function Theory.

10. NASA SP-255, Third conference on sonic boom research, NASA HEADUARTERS. —

Washington, D.C., October 29 — 30, 1970.

11. Близнюк В., Васильев Л., Вуль В., Климов В., Миронов А., Туполев А., Попов Ю., Пухов А., Черемухин Г. Правда о сверхзвуковых пассажирских самолетах. — М.: Моск. рабочий, 2000.

12. ЗавершневЮ. А., РодновА. В. Летные испытания сверхзвуковых пассажирских самолетов первого поколения по звуковому удару / Международная научно-техническая конференция «Новые рубежи авиационной науки» ASTEC'07. — Москва, 19 — 22 августа 2007.

13. Woo dward F. A. An improved method for the aerodynamic analysis of wing-body-tail configurations in subsonic and supersonic flow // NASA CR 2228, P. I, II, 1973.

Рукопись поступила 9/IX 2008 г.