Прикладная математика и механика, 2023, T. 87, № 5, стр. 869-882

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

В настоящее время одним из наиболее часто употребляемых подходов является проведение численного моделирования в ближнем возмущенном поле в рамках уравнений Эйлера или Рейнольдса, а затем использование полученных распределений избыточного давления для моделирования распространения волны звукового удара в атмосфере на основе классической квазилинейной теории и подходов геометрической акустики [4, 5] или с использованием дополненного уравнения Бюргерса [6]. Такую процедуру иногда называют двухслойной. Для корректного расчета распространения звукового удара начальное распределение следует брать на таком удалении от обтекаемого тела, где азимутальные возмущения малы относительно продольных [3, 7], т.е. начальное распределение должно быть локально осесимметричным с пренебрежимо малым поперечным течением.

Для реальных компоновок самолетов с крылом, воздухозаборниками и другими особенностями [6], способными порождать поперечное течение на значительных удалениях от тела, эти условия могут быть обычно удовлетворены на расстояниях больше пяти длин самолета L. Однако, параметрические расчеты полной трехмерной компоновки в рамках уравнений Рейнольдса с определением возмущений давления, которые малы по сравнению с величинами в набегающем потоке, на расстоянии 5–10L также является процедурой, требующей большой точности и очень существенных вычислительных ресурсов. Это связано с тем, что при удалении от летательного аппарата диссипативные свойства численной схемы приводят к увеличению погрешности расчета. Поэтому на практике часто берут избыточное давление на относительно близком расстоянии порядка 1–3L, что ведет к погрешностям в определении сигнатур звукового удара около земли из-за неучета существенного поперечного течения.

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

2. Метод мультипольного разложения. Метод мультипольного разложения [9] основан на точном решении линеаризованного уравнения для потенциала возмущенной скорости φ от сверхзвукового пассажирского самолета:

где в цилиндрической системе координат x – продольная координата, r – радиальная координата, $\theta $ – азимутальная; $\beta = \sqrt {{{M}^{2}} - 1} $, М – число Маха набегающего потока. Следует отметить, что данное уравнение предполагает малость возмущений скорости относительно скорости набегающего потока, и в нем оставлено только линейное приближение.

Решение уравнения (2.1) можно записать, разместив так называемые мультиполи с плотностью ${{A}_{n}}\left( x \right)$, ${{B}_{n}}\left( x \right)$ непрерывным образом вдоль оси х [9]:

где $\xi $ – переменная интегрирования, n – номер мультиполя, а функция интенсивности мультиполя ${{g}_{n}}$:

В характеристических переменных ($\eta = x - \beta r$)

F-функцию Уитема [3] можно определить из линеаризованного соотношения для малых возмущений давления $\Delta p$: где p_∞ – давление окружающего воздуха, γ – отношение удельных теплоемкостей. Тогда где u_∞ – скорость набегающего потока.

Следует отметить, что похожие формулы указаны, например, в [10, 11], однако там дифференцирование и интегрирование по $\eta $ проводилось иначе, что ведет к некоторым различиям.

Пусть из CFD-расчета известно поле возмущений давления $\Delta p\left( {x,R,\theta } \right)$ на некоторой цилиндрической поверхности радиуса R вокруг летательного аппарата, на которой течение еще не является локально осесимметричным, но возмущения можно считать малыми. Из соотношения (2.4) можно найти F-функцию Уитема на цилиндре R. Раскладывая ее в ряд Фурье

и вычислив интегралы, отвечающие за симметричную и несимметричную части, можно получить интегральные уравнения для ${{A}_{n}}$, ${{B}_{n}}$:

Несколько отличающиеся интегральные уравнения указаны в [10, 11], где дальнейшее их решение основано на разложении функции ${{G}_{n}}$ в ряд Тэйлора

Такое разложение становится некорректным по мере уменьшения радиуса цилиндра R, т.к. при приближении к поверхности летательного аппарата ряд начинает расходиться. Кроме того, чем больше порядок мультиполя, тем больше требуется членов ряда Тэйлора. В результате расчетов авторы [11] получили расходимость решения при уменьшении R, по мере увеличения n и J.

В настоящей работе интегральные уравнения (2.8) решаются прямым численно-аналитическим методом второго порядка точности без дополнительных упрощений и асимптотических разложений, что позволяет избежать расходимости метода [10, 11]. Для этого в интегральных уравнениях относительно неизвестной функции $A(\xi )$ (2.8) вида

производится переинтерполяция известных функций на равномерное разбиение $[{{\eta }_{1}}, \ldots ,{{\eta }_{K}}]$ с шагом $\Delta \eta $ (обычно производится разбиение не менее, чем на 500 отрезков). При этом функции предполагаются кусочно-линейными на каждом из отрезков. Тогда в каждой из точек ${{\eta }_{k}}$, $k = 2,...,K$: где ${{A}_{z}} = A({{\eta }_{z}})$, ${{H}_{k}} = H({{\eta }_{k}})$, ${{G}_{{k - z + 1}}} = G({{\eta }_{{k - z + 1}}}) = G({{\eta }_{k}} - {{\eta }_{z}})$. Обозначим где

Интегралы ${{M}_{{k,z}}}$ вычисляются аналитически:

Интегралы ${{D}_{{k,z}}}$ вычисляются аналогично. В результате получается система алгебраических уравнений, которая решается стандартным методом.

Найдя ${{A}_{n}},{{B}_{n}}$ можно получить асимптотику при $r \to \infty $, которая соответствует локально осесимметричному полю и нужна для решения задачи распространения звукового удара:

Следует отметить, что возмущение давления, соответствующее этой F-функции, приписывается обратно цилиндру радиуса R. Это связано с тем, что возмущения в ближнем поле ищутся обычно при однородном набегающем потоке, и асимптотика соответствует именно такой постановке. Но в реальности при удалении от самолета играет роль изменение атмосферы с высотой и нелинейность, которые не учитываются данной асимптотикой.

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

1) из соотношения (2.4) определить $F\left( {\eta ,\theta ,R} \right)$ на цилиндре радиуса R по результатам численного моделирования в рамках уравнений Эйлера или Рейнольдса;

2) разложить F в ряд Фурье на цилиндре по формуле (2.6);

3) вычислить интегралы (2.7);

4) найти ${{A}_{n}},{{B}_{n}}$ из интегральных уравнений (2.8);

5) определить асимптотику ${{F}^{\infty }}\left( {\eta ,\theta } \right)$ в дальнем поле по формуле (2.9);

6) найти скорректированное значение избыточного давления по формуле (2.4) и приписать его цилиндру радиуса R.

Далее скорректированную сигнатуру избыточного давления можно использовать в программах распространения звукового удара до земли, т.к. она будет соответствовать асимптотике дальнего поля, где течение локально осесимметрично (отсутствует азимутальное перетекание) и будет выполняться закон затухания $\sim {\kern 1pt} {{r}^{{ - 1/2}}}$.

По данным результатам была разработана программа MPDtool для коррекции эпюр избыточного давления с помощью метода мультипольного разложения. Ниже представлены результаты ее использования для расчета звукового удара простых тел, которые рассматривались в First Sonic Boom Prediction Workshop (SBPW) в качестве тестовых случаев [12]: SEEB-ALR – простое осесимметричное тело и DWB-69 – компоновка фюзеляж–треугольное крыло с углом стреловидности 69°.

3. Расчетные исследования обтекания осесимметричного тела. Проведены исследования звукового удара от простого осесимметричного тела SEEB-ALR [12] (рис. 1), которое моделирует носовую часть СПС. Сначала предполагалось, что модель длиной L = = 0.4491 м расположена в аэродинамической трубе [12] с нулевым углом атаки при параметрах набегающего потока М = 1.6, T_∞ = 198.43 K (температура набегающего потока), p_∞ = 16148.93 Па, что соответствует числу Рейнольдса по длине модели $\operatorname{Re} = 4.36 \times {{10}^{6}}$. Данные параметры выбраны для валидации ближнего поля возмущений модели, которое было измерено в аэродинамической трубе [12].

Проводился расчет ближнего поля в рамках уравнений Рейнольдса, замкнутых с помощью модели турбулентности Спалара–Аллмараса [13]. В расчетной области использовалась многоблочная структурированная сетка с семью миллионами ячеек, которая позволяет хорошо разрешить течение около поверхности тела. Размеры пристенных ячеек задавались так, чтобы Y+ первой ячейки по нормали к поверхности был меньше 1. Численная задача решалась стандартной программой с помощью неявного метода конечного объема сквозного счета. Схема имеет второй порядок точности по пространству. Численный алгоритм аналогичен используемому в [14].

Сравнение результатов численного моделирования в ближнем поле на расстоянии R/L = 1.2 вниз от модели с результатами эксперимента [12] приведено на рис. 2. Видно, что результаты расчета возмущений в ближнем поле лежат близко к экспериментальным данным.

Для задачи распространения звукового удара были сделаны предположения [12], что данное оживальное тело имеет длину L = 74.8 м и летит на высоте H = 16 764 м с тем же числом Маха 1.6. Ближнее поле считалось аналогично, в результате чего были получены возмущения давления на цилиндрах радиуса R/L = 0.25, 0.5, 1, 1.2, 2.38.

Затем избыточное давление на этих цилиндрах корректировалось с помощью программы MPDtool на основе мультипольного разложения. Дискретизация по азимутальному углу составляла 2°. После этого полученная эпюра задавалась в качестве исходных данных для задачи распространения звукового удара до земли, которая решалась в рамках классического подхода [4, 5] с помощью собственной программы sbProp.

Результаты расчетов сигнатуры звукового удара на земле от осесимметричного тела SEEB-ALR представлены на рис. 3 для наибольшего и наименьшего из рассмотренных цилиндров R/L = 2.38 и 0.25. Видно, что если не использовать мультипольное разложение, то эпюры на земле заметно отличаются для случаев расчета от ближнего поля на R/L = 2.38 (кривая 3 на рис. 3) и R/L = 0.25 (кривая 1 на рис. 3). Кривая 1 лежит заметно ниже под плато и отличается перед задним скачком. С увеличением R/L результаты сходятся к корректному решению, которое не зависит от R/L.

Использование мультипольной коррекции с R/L = 2.38 (кривая 4 на рис. 3) ведет к практически такой же сигнатуре на земле, как и без мультипольной коррекции (кривая 3). Это говорит о том, что кривая 3 практически совпадает с корректным (предельным) решением, т.к. толщина рассматриваемого тела мала, крыла нет и возмущения осесимметричны. Поэтому сигнатура с R/L = 2.38 без мультипольной коррекции и с ней практически сливаются.

Применение мультипольной коррекции для эпюры избыточного давления с R/L = = 0.25 и последующее решение задачи распространения звукового удара до земли приводят к кривой 2, которая ложится на кривую 4 с R/L = 2.38. Это говорит о возможности считать ближнее поле возмущений от этого тела до R/L = 0.25 и применять его для задачи распространения звукового удара.

Следует отметить, что данное тело является практически осесимметричным (с точностью до изготовления [12]). Поэтому ближнее поле осесимметрично и мультипольное разложение по азимутальному углу не дает коррекции. Однако коррекция происходит из-за разложения по r, т.к. ближнее поле еще не соответствует асимптотике $\sim {\kern 1pt} {{r}^{{ - 1/2}}}$, которая используется при решении задачи распространения звукового удара до земли.

4. Расчетные исследования компоновки крыло–фюзеляж. Проведены исследования звукового удара от компоновки DWB-69 фюзеляж–треугольное крыло с углом стреловидности 69° [12] (рис. 4). Крыло имеет ромбовидную форму с толщиной 5%. Сначала предполагалось, что модель длиной L = 0.1752 м расположена в аэродинамической трубе [12] с нулевым углом атаки при параметрах набегающего потока М = 1.7, T₀ = = 308.2 K (полная температура набегающего потока), p₀ = 110235 Па (полное давление). Как и в случае оживального тела, проведена валидация ближнего поля модели [12], которое рассчитывалось для половины симметричной модели аналогично предыдущему пункту. Использовалась многоблочная структурированная сетка с пятью-десятью миллионами ячеек с Y+ первой ячейки по нормали к поверхности меньше 1.

На рис. 5 представлено сравнение численного решения в рамках уравнений Рейнольдса (кривая 3) на расстоянии R/L = 3.6 вниз от модели с результатами эксперимента [12] (сплошная кривая 1 и разброс экспериментальных данных, обозначенный серой областью 2). Результаты расчета (кривая 3) ближнего поля лежат близко к экспериментальным данным и попадают в разброс расчетных значений других авторов [12].

Для задачи распространения звукового удара были сделаны предположения [12], что данная компоновка имеет длину L = 26.95 м и летит на высоте H = 16764 м с тем же числом Маха 1.7. Ближнее поле считалось аналогично, в результате чего были получены возмущения давления на цилиндрах радиуса R/L = 0.5, 1, 2, 3, 3.6.

Затем, как и ранее, избыточное давление на этих цилиндрах корректировалось с помощью программы MPDtool на основе мультипольного разложения, а полученная эпюра задавалась в качестве исходных данных для задачи распространения звукового удара до земли по программе sbProp. Количество мультиполей в разложении было n = 10.

Результаты расчетов сигнатуры звукового удара на земле от компоновки DWB-69 представлены на рис. 6 для расчетов от R/L = 0.5, 1, 2, 3.6 без использования мультипольной коррекции. Эпюры сходятся к предельному корректному значению при увеличении R/L, как это хорошо видно на промежуточном и замыкающем скачке (передний скачок для всех эпюр совмещен с началом координат).

Использование мультипольной коррекции (рис. 7) ведет к тому, что на земле все сигнатуры (при расчете от R/L = 0.5, 1, 2, 3.6) практически сливаются. Кроме того, эпюра от R/L = 3.6 без коррекции (кривая 4 на рис. 6) очень близка к эпюрам после коррекции (все кривые на рис. 7). Это говорит о сходимости метода и возможности использования ближнего поля от R/L = 0.5 с применением мультипольного разложения.

5. Применение интеграла Жилина. Ранее в [15] был разработан подход для определения звукового удара на земле от СПС с использованием интеграла Жилина. В этом случае ближнее поле может быть рассчитано в рамках, например, уравнений Эйлера или Рейнольдса. Предполагается, что возмущения, вносимые летательным аппаратом в набегающий поток, малы, и используется линейное приближение.

В [15] проведена контрольная поверхность из плоскостей S₁–S₄ (рис. 8). S₁ отделяет возмущенный поток от невозмущенного (конус Маха). S₂ обычно параллельна плоскости XZ, лежит ниже тела и не пересекает его. Точки А и В лежат на линии пересечения плоскости S₂ и поверхности S₁, S₃ и S₄ – огибающие обратных конусов Маха, проходящие через отрезок АВ.

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

где $u$, ${v}$ – возмущения продольной и вертикальной компоненты скорости соответственно; x₀ – координата точки пересечения S₃ с осью x. Выписанные формулы верны для направления под самолетом, а для других азимутальных углов они усложняются. Также можно записать: где интегралы вычисляются вдоль отрезка АВ.

Далее для расчетов будет использоваться второй из этих интегралов относительно ${v}$. Тогда F-функцию Уитема можно вычислить по стандартной формуле [3]:

Для сравнения такого подхода с предыдущими проведены расчеты звукового удара от компоновки DWB-69 с длиной L = 26.95 м, которая летит на высоте H = 16 764 м с числом Маха 1.7 (как в разд. 3). Расчет ближнего поля представлен в разд. 3. В результате были получены возмущения давления в плоскостях H/L = 0.3, 0.5, 1, 2, 3. Затем при помощи второго интеграла (5.1) находилась F-функция Уитема, которая использовалась для решения задачи распространения звукового удара до земли.

Результаты расчетов сигнатуры звукового удара на земле представлены на рис. 9 для разных вариантов H/L. Эпюры сходятся к предельному значению при увеличении H/L, как это видно на промежуточном скачке (передний скачок для всех эпюр совмещен с началом координат). Однако, избыточное давление для H/L = 0.3 заметно отличается от случая H/L = 3.

Таким образом, в [15] были получены интегральные соотношения для возмущенных скоростей, которые можно использовать для задачи распространения звукового удара до земли. Но такой подход, в целом, имеет схожие недостатки, что и двухслойный подход, указанный во введении. Как указано в [15], способом решения задачи может служить прямой расчет пространственного обтекания тела с выходом на асимптотическое решение, где оно развивается в каждой меридиональной плоскости независимо с зависимостью $\sim {\kern 1pt} {{r}^{{ - 1/2}}}$. Такое поведение для сложных компоновок достигается обычно после 5 длин летательного аппарата.

Другой вариант, указанный, например, в [3] содержит расчет ближнего поля в рамках уравнений Эйлера с определением интеграла Жилина на плоскости H/L = 0.3–0.5, затем расчет распространения избыточного давления по закону Келлера [16] до ~5L. После чего проводится расчет распространения звукового удара до земли в атмосфере по классической теории [4].

Таким образом, для надежного предсказания эпюры избыточного давления на земле исходные данные для задачи распространения звукового удара надо брать на расстоянии больше 5L даже в случае использования интеграла Жилина.

Заключение. Разработан метод расчета звукового удара от СПС с использованием мультипольного разложения с численно-аналитическим прямым способом решения возникающих интегральных уравнений второго порядка точности. Метод не имеет расходимости даже на малых расстояниях от летательного аппарата ~0.25 его длины.

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

Авторы благодарят Ю.С. Суворову и К.Г. Хайруллина за помощь при проведении расчетов.

Публикация подготовлена в рамках реализации Программы создания и развития научного центра мирового уровня “Сверхзвук” на 2020–2025 годы при финансовой поддержке Министерством науки и высшего образования Российской Федерации (соглашение о предоставлении гранта в форме субсидий из федерального бюджета на осуществление государственной поддержки создания и развития научных центров мирового уровня, выполняющих исследования и разработки по приоритетам научно-технологического развития от 17 мая 2022 г. № 075-15-2022-1023).