Журнал вычислительной математики и математической физики, 2021, T. 61, № 11, стр. 1904-1926
Постановка и метод решения задачи бимформинга для локализации акустического источника на основе данных вычислительного эксперимента
Т. К. Козубская 1, *, Г. М. Плаксин 1, **, И. Л. Софронов 2, 3, ***
1 ИПМ РАН
125047 Москва, Миусская пл., 4, Россия
2 МФТИ
141700 М.o., Долгопрудный, Институтский пер., 9, Россия
3 “Шлюмберже”
125171 Москва, Ленинградское ш., 16а, стр. 3, Россия
* E-mail: tatiana.kozubskaya@gmail.com
** E-mail: glebyp@gmail.com
*** E-mail: ilsofronov@gmail.com
Поступила в редакцию 11.06.2020
После доработки 11.06.2020
Принята к публикации 07.07.2021
Аннотация
Статья посвящена постановке задачи бимформинга, целью которой является определение непрерывной функции акустического источника на основе анализа пространственно-временных данных вычислительного эксперимента в задачах аэродинамики и аэроакустики, и разработке численного метода ее решения. Под вычислительным экспериментом понимается максимально точное воспроизведение исследуемого турбулентного течения путем его численного моделирования, что сегодня становится возможным благодаря использованию вихреразрешающих подходов, методов повышенной точности, сеток большой размерности и мощной суперкомпьютерной техники. В работе сформулированы условия на параметры сеток дискретизации области источника монопольного типа и сетки микрофонов для обеспечения корректности и точности работы алгоритма. Предложен способ задания входящих в указанные условия численных параметров, корректность которых контролируется величиной числа обусловленности оператора бимформинга. В случае нарушения условий корректности рассмотрена регуляризация оператора бимформинга с целью получения функции источника с приемлемой точностью. Показано, что предложенный численный метод бимформинга с высокой точностью восстанавливает функцию источника на тестовых примерах, а его применение к расчетным данным вычислительного эксперимента по турбулентному обтеканию сегмента крыла с механизацией дает хорошее согласие с полученными в ходе эксперимента результатами анализа акустического поля. Библ. 20. Фиг. 23. Табл. 5.
1. ВВЕДЕНИЕ
В аэроакустике под бимформингом понимают технологию нахождения акустического источника для волнового уравнения, в предположении, что он приближенно воссоздает звуковое поле, формируемое летящим (или обдуваемым) объектом. Целью бимформинга является определение мощности и спектрального состава этого источника в анализируемой области пространства на основе данных, полученных на расположенных вдали микрофонах. В современных экспериментальных исследованиях бимформинг представляет собой хорошо разработанную технологию, предоставляющую информацию высокой практической значимости с точки зрения разработки малошумных конструкций летательных аппаратов.
Потребность в выявлении источников чрезмерного шума в технических изделиях (летательные аппараты, автомобили, поезда, строительная техника и т.д.) по итогам акустических измерений постоянно растет. Для решения этой задачи наибольшее распространение получили подходы, основанные на использовании так называемой фазированной микрофонной решетки (phased microphone arrays) [1], [2], которые используются для локализации источника шума в натурных и лабораторных экспериментах. В основе таких подходов лежит когерентное суммирование значений на микрофонах, чтобы, с одной стороны, усилить акустический сигнал, излучаемый из точки фокуса микрофонной решетки, и, с другой стороны, ослабить сигнал, излучаемый из точек вне фокуса. Выполнять суммирование можно различными способами, например, с задержкой, вызванной разницей времени прихода сигнала от точки фокуса к точкам расположения микрофонов. Применение данных методов подразумевает наличие априорной информации о ширине спектра и характере излучения, а также о когерентности и компактности источников.
Для повышения точности решения задачи бимформинга, особенно на низких частотах, используются адаптивные алгоритмы бимформинга [3], [4], однако, эти методы на практике оказываются неустойчивыми к возмущениям акустических измерений, что затрудняет аккуратную оценку интенсивности источников. Та же проблема возникает и при использовании метода CLEAN [5]. Для устранения этого недостатка были предложены методы DAMAS [6] и DAMAS2 [7], основанные на уточнении карт шума, полученных методом conventional beamforming [1] за счет выделения нескоррелированного распределения источников путем итерационной деконволюции с так называемой point spread function, описывающей излучение монополя. Были также предложены модификации DAMAS-C [8] и CLEAN-SC [9] для решения задачи в случае наличия когерентных источников шума. Ограничением данных методов является использование монопольных излучателей, в то время как во многих случаях источниками шума являются мультиполи.
В [10] предлагается подход, альтернативный вышеописанным. Он предназначен как для нескоррелированных, так и скоррелированных источников шума. В основе метода лежит псевдообращение прямоугольной матрицы переноса излучения от источников к микрофонам. Задача бимформинга решается итерационно; на каждом шаге удаляются источники с малой интенсивностью, что приводит к переопределенности задачи, начиная с некоторого этапа. Однако данный подход не дает возможности восстанавливать гладкую функцию источника и может вызвать трудности при использовании данных, полученных в ходе численного эксперимента в аэроакустических задачах.
Мы упомянули лишь очень малую часть работ по методам бимформинга; более подробную информацию можно найти, например, в обзорной статье [11]. Практически все рассмотренные в ней методы ориентированы на экспериментальные исследования с присущими им преимуществами и ограничениями. С ростом производительности современных суперкомпьютеров и развитием численных алгоритмов резко выросла роль вычислительного эксперимента как еще одного эффективного средства инженерных исследований. При этом алгоритм бимформинга применительно к вычислительному эксперименту модифицируется, приобретает свои характерные особенности и, в том числе, свои преимущества. Преимущества, главным образом, связаны с возможностью использования всей объемной информации, предоставляемой вычислительным экспериментом.
В предпринятом исследовании мы формулируем математическую задачу нахождения функции акустического источника с носителем в предписанной области пространства, опираясь на возможность обработки большого объема пространственно-временных данных, получаемых в ходе численного моделирования. При этом мы рассматриваем задачи внешнего обтекания, численно решаемые с помощью вихреразрешающих подходов, к которым относятся: прямое численное моделирование (DNS, Direct Numerical Simulation), моделирование крупных вихрей (LES, Large Eddy Simulation) и гибридные методы, сочетающие подход LES с решением осредненных по Рейнольдсу уравнений Навье-Стокса (RANS, Reynolds Averaged Navier-Stokes) в пристеночных областях. Задача бимформинга использует данные об акустическом давлении на некоторой поверхности с виртуально заданными микрофонами, охватывающей область предполагаемого акустического источника. Анализ формулируемой в разд. 2 задачи проводится в частотной области, переход в которую осуществляется преобразованием Фурье исходных данных вычислительного эксперимента. В рассматриваемом нами подходе звуковое давление удовлетворяет уравнению Гельмгольца в движущейся среде в трехмерном пространстве с непрерывной правой частью, соответствующей заранее неизвестной функции искомого источника в ограниченной области на некоторой выбранной поверхности. В разд. 3 мы проводим дискретизацию задачи на основе конечно-элементного подхода и интегрального представления решения, используя кусочно-линейный базис для функции источника. Раздел 4 статьи посвящен выводу физически обоснованных условий для выбора сетки для дискретизации области источника и точек расположения микрофонов (иначе говоря, сетки микрофонов) с целью построения корректной матрицы бимформинга, обеспечивающей устойчивое решение с максимально возможной точностью. Показано, что параметры сеток источника и микрофонов всецело определяются длиной волны. В разд. 5 мы описываем приемы получения конкретных значений для шагов сеток источника и микрофонов на основе специально подобранных тестовых задач. В результате определяются достаточно универсальные интервалы значений отношения шагов сетки к длине волны для процедуры выбора оптимальных параметров метода. Оценивается также точность решения на примере тестовых функций. Кроме того, рассматривается и обосновывается подход применения регуляризации по Тихонову задачи бимформинга для случаев, когда нарушаются условия построения корректной матрицы. Раздел 6 содержит результаты расчетов по восстановлению некоторых заранее заданных тестовых функций источника при использовании поверхности микрофонов, соответствующей геометрической конфигурации вычислительного эксперимента. В разд. 7 рассматриваются данные решения задачи обтекания крылового сегмента с выпущенной механизацией и приводятся результаты восстановления интенсивности и формы источников звука разработанным подходом для низко- и высокочастотных третьоктавных полос спектра. Получаемые решения задачи бимформинга хорошо согласуются с мгновенным распределением ближнего поля давления задачи обтекания. Отметим, что особенностью представленного в статье нового подхода является возможность получения корректной матрицы бимформинга без необходимости привлечения приемов регуляризации при формулировке и (или) решении задачи.
2. НЕПРЕРЫВНАЯ ЗАДАЧА БИМФОРМИНГА
2.1. Постановка задачи
В линейном приближении аэроакустики [14], лежащем в основе бимформинга, звуковое давление$~$ $p$ удовлетворяет задаче Коши для волнового уравнения в неподвижной среде:
(1)
$\begin{gathered} \frac{1}{{{{c}^{2}}}}\frac{{{{\partial }^{2}}p}}{{\partial {{t}^{2}}}} + 2~~\frac{{{{\operatorname{M} }_{\infty }}}}{c}\frac{{{{\partial }^{2}}p}}{{\partial t\partial {{x}_{1}}}} - \Delta p + ~\operatorname{M} _{\infty }^{2}\frac{{{{\partial }^{2}}p}}{{\partial x_{1}^{2}}} = q(x,~t),\quad ~x \in {{\operatorname{R} }^{3}},\quad t > 0, \\ p(x,~0) = \frac{{\partial p}}{{\partial t}}(x,0) = 0,\quad x \in ~{{\operatorname{R} }^{3}}{\kern 1pt} . \\ \end{gathered} $Фиг. 1.
Поверхность сегмента крыла и поверхность распределенного акустического источника $S$ в двумерном (а) и трехмерном (б) случаях.
Пусть существует поверхность D, не имеющая пересечений с поверхностью S, $S \cap D = \emptyset $, и заданы значения $p(x,t)$, $x \in D$, на этой поверхности. Поверхность D соответствует предполагаемому расположению микрофонов в задаче бимформинга.
Тогда задача бимформинга ставится следующим образом: по значениям ${{\left. p \right|}_{D}}$ акустического поля, удовлетворяющего уравнению (1), требуется восстановить функцию источника $q(x,t)$.
Будем рассматривать поставленную задачу в случае гармонических колебаний по времени:
где f – частота звука. Подставив эти функции в (1), получим уравнение в частотной области для функции $P(x)$ в движущейся среде:(2)
${{k}^{2}}P + {{\Delta }}P - 2ik{{\operatorname{M} }_{\infty }}\frac{{\partial P}}{{\partial {{x}_{1}}}} - \operatorname{M} _{\infty }^{2}\frac{{{{\partial }^{2}}P}}{{\partial x_{1}^{2}}} = - Q(x),\quad k = \frac{{2\pi f}}{c}.$Таким образом, будем решать следующую задачу бимформинга в частотной области: по значениям ${{\left. P \right|}_{D}}$ на поверхности D акустического поля, удовлетворяющего уравнению (2), необходимо восстановить значения функции $Q(x)$ на поверхности $S$.
Замечание. В вычислительной аэроакустике исследуется набор таких частотных задач для различных значений параметра $f$; переход в частотную область осуществляется путем применения (быстрого) преобразования Фурье к временным массивам давления, полученных в расчете в точках расположения виртуальных микрофонов.
2.2. Интегральное представление звукового поля
Известно, что уравнение (2) имеет следующее решение [15]:
(3)
$\begin{gathered} P\left( x \right) = \mathop \smallint \limits_S \,Q(y){{G}_{{{{\operatorname{M} }_{\infty }}}}}(x - y)dy,\quad x \in D, \\ {{G}_{{{{{\text{M}}}_{\infty }}}}}(x) = \frac{1}{{4\pi }}\frac{{{{e}^{{ - ik{\kern 1pt} {\text{'}}(x{\kern 1pt} {\text{'}} - {{\operatorname{M} }_{\infty }}{{x}_{1}})}}}}}{{x{\kern 1pt} {\text{'}}}}, \\ x{\text{'}} = \sqrt {x_{1}^{2} + (1 - \operatorname{M} _{\infty }^{2})(x_{2}^{2} + x_{3}^{2})} ,\quad k{\text{'}} = \frac{k}{{1 - \operatorname{M} _{\infty }^{2}}}. \\ \end{gathered} $3. ДИСКРЕТНАЯ ЗАДАЧА БИМФОРМИНГА
3.1. Постановка задачи
Пусть выполнены следующие условия:
1) на поверхности S задано N точек $\{ {{y}^{n}}\} _{{n = 1}}^{N}$, которые формируют сетку источника;
2) на поверхности D задано M точек $\{ {{x}^{m}}\} _{{m = 1}}^{M}$, которые назовем микрофонами, причем количество точек сетки источника $N$ не превышает количество микрофонов M.
Введем следующие обозначения:
• ${\mathbf{s}} = {{({{s}_{1}},{{s}_{2}},...,{{s}_{N}})}^{{\text{т}}}},\;{{s}_{n}} = Q({{y}^{n}})$, – дискретный вектор интенсивности источника,
• ${\mathbf{d}} = {{({{d}_{1}},{{d}_{2}},...,{{d}_{M}})}^{{\text{т}}}},\;{{d}_{m}} = P({{x}^{m}})$, – вектор звуковых сигналов, полученных микрофонами,
• ${{\mathcal{T}}_{a}}$ – некоторая дискретная аппроксимация оператора $\mathcal{T}$ из (4) (см. п. 3.2).
Сформулируем задачу нахождения интенсивности распределенного звукового источника следующим образом: необходимо найти такой вектор ${\mathbf{s}}$, что выполняется соотношение
(5)
${\mathbf{d}} = {{\mathcal{T}}_{a}}{\mathbf{s}},\quad {{\mathcal{T}}_{a}} \in {{\mathbb{R}}^{{M \times N}}}.$(6)
$\left\| {{\mathbf{\tilde {d}}} - {{\mathcal{T}}_{a}}{\mathbf{s}}{\kern 1pt} } \right\|_{2}^{2} \to \mathop {{\text{min}}}\limits_{\mathbf{s}} ,$3.2. Метод решения задачи бимформинга
Известно, что при $M > N$ задача (6) имеет решение (которое можно найти методом наименьших квадратов):
Введем на поверхности S сетку с узлами в $\{ {{y}_{n}}\} _{{n = 1}}^{N}$ для кусочно-линейного восполнения функций, непрерывных на S. Пусть $\{ {{\psi }_{n}}(y)\} _{{n = 1}}^{N}$ – соответствующие базисные элементы пространства кусочно-линейных функций. Тогда вместо функции источника $Q(y)$ будем искать ее кусочно-линейную аппроксимацию
Подставив эту аппроксимацию в (3), получим следующие выражения для звуковых сигналов на поверхности D:(7)
${{\mathcal{T}}_{{mn}}} = \mathop \smallint \limits_S \,{{\psi }_{n}}(y){{G}_{{{{\operatorname{M} }_{\infty }}}}}({{x}_{m}} - y)dy.$Рассмотрим подробнее конструкцию базисных функций и способ вычисления интеграла (7).
3.2.1. Расположение источника на линии. Иногда оказывается достаточным располагать источник на линии, а не на поверхности. Это бывает, например, когда проводится бимформинг для фрагментов конструкций на частотах с длиной волны порядка одного из размеров поверхности S. Пример такой задачи приведен в п. 7.1, где рассматривается случай обтекания крылового сегмента достаточно малого поперечного размера. Тогда становится возможным аппроксимировать поверхность источника $S$ некоторой ломаной линией, на которой заданы узлы $\{ {{X}^{n}}\} _{{n = 1}}^{N}$ одномерной сетки источника. Базисные функции ${{\psi }_{n}}$ для кусочно-линейного восполнения имеют простейшую треугольную форму (см. фиг. 2):
3.2.2. Расположение источника на поверхности. В общем случае для аппроксимации поверхности источника S будем использовать треугольную сетку и соответствующую кусочно-треугольную поверхность. Узлы сетки источника при этом размещаются в вершинах треугольников, а базисные функции ${{\psi }_{n}}$ для двумерного кусочно-линейного восполнения имеют пирамидальную форму.
3.2.3. Вычисление элементов матрицы бимформинга ${{\mathcal{T}}_{a}}.~$ Аппроксимация интегралов (7) осуществляется по формулам Гаусса: квадратурами на отрезках $[{{X}^{{n - 1}}},{{X}^{n}}]$ для одномерных базисных функций и кубатурами на треугольниках для двумерных базисных функций (см., например, [13]). Количество узлов k в носителе базисной функции, определяющее точность вычисления ${{\mathcal{T}}_{{mn}}}$, существенно влияет на получаемое решение задачи бимформинга. Оптимальной по соотношению точности и вычислительных затрат является формула с k = 4 для отрезка, как показывают результаты теста, приведенного ниже в п. 6.1, и k = 7 для треугольника в аналогичном двумерном тесте (см. п. 6.2).
4. ПОСТРОЕНИЕ КОРРЕКТНОЙ МАТРИЦЫ БИМФОРМИНГА
Формулировка дискретной задачи бимформинга содержит два множества точек – сетку для функции источника и сетку виртуальных микрофонов. Параметры этих сеток влияют на число обусловленности ${{C}_{a}}~$ матрицы $\mathcal{T}_{a}^{{\text{*}}}{{\mathcal{T}}_{a}}$ и на погрешность ε вычисления функции источника. Численные эксперименты показывают, что даже незначительное изменение в величине шага сетки для источника может на порядки изменить величину ${{C}_{a}}~$ (см. разд. 5). Поэтому мы будем говорить, что матрица построена корректно, если ${{C}_{a}}~$ не является чрезмерно большим (например, не превышает 105), а ε удовлетворяет желаемой величине на характерных тестовых функциях источника. Для краткости будем называть такие матрицы корректными.
Здесь мы приводим некоторые физические соображения, позволяющие получить оценки параметров сеток источника и микрофонов, чтобы матрица была корректной (хотя данные рассуждения не являются строгими, они дают понимание о порядке значений нужных нам величин, участвующих в формировании матрицы ${{\mathcal{T}}_{a}}$). Затем, в разд. 5, мы предлагаем алгоритм подбора конкретных значений сеточных параметров для окончательной формулировки дискретной задачи с заданной геометрией поверхностей S функции источника и D микрофонов.
4.1. Сетка для аппроксимации функции источника
Прежде всего, потребуем, чтобы соседние точки сетки были достаточно различимы с точки зрения вклада значений в них функции источника в акустическое поле. Это означает, что сеточные узлы не должны располагаться на расстоянии, существенно меньшем длины звуковой волны λ для рассматриваемой частоты. Данное ограничение связано с дифракционным пределом и может быть проиллюстрировано примером, изображенным на фиг. 3, где distS – шаг сетки функции источника. Очевидно, что должно выполняться условие ${{\operatorname{dist} }_{S}} \gtrsim \lambda $. С другой стороны, для обеспечения как можно меньших значений ε нельзя располагать узлы сетки с шагом, существенно превышающим длину волны, т.е. необходимо выполнение условия ${\text{dis}}{{{\text{t}}}_{S}} \lesssim \lambda $. Таким образом, мы приходим к первому условию: шаг сетки представления функции источника должен быть сопоставим с длиной волны:
4.2. Сетка микрофонов
Поверхность микрофонов D характеризуется двумя параметрами: расстоянием distSD от поверхности S и расстоянием между соседними микрофонами distD. На фиг. 4 показаны две соседние точки сетки s1 и s2 на S, ${\text{dis}}{{{\text{t}}}_{S}} = \left| {{{s}_{1}}{{s}_{2}}} \right|$, два ближайших к ним микрофона m1 и m2, шаг сетки микрофонов ${\text{dis}}{{{\text{t}}}_{D}} = \left| {{{m}_{1}}{{m}_{2}}} \right|$. Для упрощения проведения оценок считаем, что S и D параллельные плоскости, и ${\text{dis}}{{{\text{t}}}_{{SD}}} \gg {\text{dis}}{{{\text{t}}}_{S}}$.
Фиг. 4.
Схема расположения соседних узлов сетки источника s1 и s2 на поверхности S и микрофонов m1 и m2 на поверхности D.
Каждый из микрофонов m1 или m2 принимает сигнал от функции источника на отрезке [s1, s2]. Для различения значений функции в этих точках сетки необходимо, чтобы расстояния от какого-либо из микрофонов до каждой из них заметно различались по отношению к длине волны, т.е. разность расстояний составляла величину const λ независимо от distSD. Для изображенной на фиг. 4 конфигурации видно, что для первого микрофона $\left| {{{s}_{1}}{{m}_{1}}} \right| \approx \left| {{{s}_{2}}{{m}_{1}}} \right|$, поэтому потребуем выполнение указанного условия для второго микрофона:
(9)
$\left| {{{s}_{1}}{{m}_{2}}} \right| - \left| {{{s}_{2}}{{m}_{2}}} \right| = \operatorname{const} \lambda .$(10)
${\text{dis}}{{{\text{t}}}_{D}} \approx {\text{dis}}{{{\text{t}}}_{{SD~}}}\sin \phi \approx {\text{dis}}{{{\text{t}}}_{{SD}}}\frac{{\operatorname{const} \lambda }}{{{\text{dis}}{{{\text{t}}}_{S}}}},$Для оценки расстояния distSD заметим, что необходимо, чтобы микрофоны хорошо разрешали функцию на [s1, s2] даже в том неблагоприятном случае, когда в этих двух точках источник имеет одинаковую амплитуду, но противоположный знак. С учетом (8) моделью данной ситуации может служить диполь с базой ~λ. Несложные вычисления для фундаментального решения (3) при ${\text{dis}}{{{\text{t}}}_{{SD}}} \gg {\text{dis}}{{{\text{t}}}_{S}}$ показывают, что сигнал на микрофоне m1 от такого диполя с единичной амплитудой будет ${{d}_{1}}\sim {{\left( {\frac{\lambda }{{{\text{dis}}{{{\text{t}}}_{{SD}}}}}} \right)}^{2}}$. Это означает, что матрица $\mathcal{T}_{a}^{{\text{*}}}{{\mathcal{T}}_{a}}$ будет иметь маленькие собственные числа порядка $d_{1}^{2}$. Соответственно, рост числа обусловленности ${{C}_{a}}~$ матрицы $\mathcal{T}_{a}^{{\text{*}}}{{\mathcal{T}}_{a}}$ можно оценить как
Таким образом, для умеренного числа обусловленности необходимо, чтобы расстояние между поверхностями D и S не превышало нескольких десятков, а то и единиц, длин волн, т.е. Наконец, для устойчивой разрешимости задачи (6) нужно потребовать, чтобы на один узел сетки источника приходилось несколько микрофонов, т.е.Полученные оценки (8), (10)–(13) важно учитывать при построении геометрии и оценке вычислительных ресурсов, необходимых при планировании задачи бимформинга в вычислительном эксперименте. Можно сделать следующие выводы качественного характера, полезные для выбора геометрических параметров бимформинга:
1) расстояние между поверхностями D и S, а также шаги сеток на них зависят от частоты исследуемого акустического сигнала, а именно пропорциональны длине волны $\lambda $,
2) поверхность D должна “охватывать” S, но необязательно быть замкнутой.
Алгоритм нахождения конкретных значений шага сетки источника distS, шага сетки микрофонов distD и расстояния между поверхностями источника и микрофонов distSD рассматриваются в следующем разделе.
5. ВЫБОР ПАРАМЕТРОВ ДЛЯ БИМФОРМИНГА
Основным входным параметром метода является длина волны $\lambda = \frac{c}{f}(1 - \operatorname{M} _{\infty }^{2})$ исследуемого акустического сигнала.
Для проведения на модельных задачах количественных оценок параметров, нужных для корректного бимформинга, введем подборочные коэффициенты aSD, aS, aD, которые определяют расстояние между поверхностями функции источника и микрофонов, а также шаги сеток источника и микрофонов соответственно:
5.1. Одномерный случай: источник на линии
Сначала рассмотрим случай расположения источника и микрофонов на линиях. Для небольшого участка геометрии полагаем, что источник и микрофоны расположены на параллельных отрезках S и D, образуя симметричную конфигурацию (фиг. 4).
Задаем, например, aS = 1, N = 10, M = 15 и подыскиваем значения параметра aD в пределах 0 < < aD < 1, обеспечивающие умеренную величину числа обусловленности ${{C}_{a}}~$, например, ${{C}_{a}}~$ < 105. Для этого рассматриваем набор значений параметра aSD = 3, 5, 10, 20, 30, строим сетки микрофонов и вычисляем соответствующие матрицы ${{\mathcal{T}}_{a}}$. Типичные графики зависимости ${{C}_{a}}~$ от aD при различных значениях изображены на фиг. 5а. Из них видно, что подходящим является интервал значений 0.08 < aD < 0.3. Выберем значение aD = 0.16 из этого интервала и теперь подберем aS в пределах 0.3 < aS < 2, которые также обеспечивают умеренную величину числа ${{C}_{a}}~$. Типичные графики изображены на фиг. 5б, откуда видно, что величину aS можно взять из диапазона 0.4 < aS < 1.5.
Фиг. 5.
Источник на линии: зависимость числа обусловленности Ca от коэффициентов aD (а) и aS (б) при различных значениях aSD.
Для оценки точности восстановления источника на получаемых сетках можно рассмотреть, например, в качестве тестовой функции источника Гауссиан с полушириной 2λ:
5.2. Двумерный случай: источник на поверхности
Такой же анализ проведем для случая поверхностей S и D, полагая их квадратными областями, расположенными на параллельных плоскостях. Соответственно, рассматриваем квадратные сетки источника и микрофонов с пробными значениями N = 102, M = 152. Ситуация с поверхностями оказывается полностью аналогичной ситуации с линиями. Сначала находим подходящие интервалы для коэффициентов aD и aS при допустимых значениях числа обусловленности (см. фиг. 7): $0.07 < {{a}_{D}} < 0.2$, $0.5 < {{a}_{S}} < 1.3$. Затем, зафиксировав для примера значение ${{a}_{S}} = 0.8$, уточняем интервал для aD на графиках погрешности для тестовых функций (см. фиг. 8). Видим, что погрешность составляет несколько процентов в интервале $0.08 < {{a}_{D}} < 0.15$ (отметим, что из-за использования более грубой сетки, здесь рассматриваются гауссианы с полушириной 3λ и 4λ).
Фиг. 8.
C-норма ошибки $\varepsilon _{{2D}}^{1}~$и $\varepsilon _{{2D~}}^{2}$в зависимости от коэффициента aD, полуширина гауссиана 3λ (а) и 4λ (б).
Резюмируя результаты исследований, представленных в п. 5.1, 5.2, получаем следующие подходящие нам значения коэффициентов:
• для одномерного случая (источника на линии): $0.4 < {{a}_{S}} < 1.5$, $0.08 < {{a}_{D}} < 0.15$ (можно, например, зафиксировать ${{a}_{S}} = 0.6$, ${{a}_{D}} = 0.12$);
• для двумерного случая (источник на поверхности): $0.6 < {{a}_{S}} < 1.5$, $0.08 < {{a}_{D}} < 0.15$ (можно, например, зафиксировать ${{a}_{S}} = 0.8$, ${{a}_{D}} = 0.1$).
При этом численные оценки точности, проведенные на примере функции источника в виде гауссиана с полушириной 3λ, дают примерно ${{\varepsilon }_{{1D}}} = 0.02{\kern 1pt} - {\kern 1pt} 0.03$, ${{\varepsilon }_{{2D}}} = 0.04{\kern 1pt} - {\kern 1pt} 0.07$.
5.3. О построении поверхностей S и D и сеток
Для конкретной геометрии, например, сегмента крыла самолета, вначале задаем поверхность (или линию) источника S, состоящую из плоских (или прямых) фрагментов (см. фиг. 1). Затем, выбрав желаемое значение ${{a}_{{SD}}} \in [3,30]$, строим поверхность (или линию) микрофонов D, которая отстоит от S на расстояние примерно aSDλ. Геометрия D может быть достаточно произвольной с точки зрения криволинейности, поскольку на ней нам нужны только точки сетки микрофонов (в отличие от поверхности S, на которой строятся базисные функции для аппроксимации функции источника).
На поверхностях (или линиях) S и D строим сетки источника и микрофонов с узлами, примерно равноотстоящими по координатным линиям на ${\text{dis}}{{{\text{t}}}_{S}} = {{a}_{S}}~\lambda $ и ${\text{dis}}{{{\text{t}}}_{D}} = \frac{{{{a}_{D}}}}{{{{a}_{S}}}}{\text{\;dis}}{{{\text{t}}}_{{SD}}}$ друг от друга соответственно. При этом важно, чтобы distS было не меньше дифракционного предела, оцениваемого как 0.6λ и 0.4λ для поверхности и линии S соответственно.
В то же время аналогичной нижней границы для distD нет, поскольку рост ошибок ${{\varepsilon }_{{2D}}}$ и ${{\varepsilon }_{{1D}}}$ при значениях коэффициента aD, выходящих за пределы левой границы отрезка из вышеприведенных оценок (0.03 и 0.07) связан, прежде всего, с фиксацией значений M = 152 и M = 15 в модельных задачах. Если при уменьшении aD увеличивать число узлов M так, чтобы поверхность D всегда “охватывала” S, то матрица ${{\mathcal{T}}_{a}}$ остается корректной, и потери точности не происходит (как показывают дополнительные численные эксперименты).
5.4. Регуляризация задачи бимформинга для некорректной матрицы
В реальных расчетах могут возникать нарушения условий (8), (10)–(13) на шаги сеток функции источника и микрофонов и их удаленности друг от друга. Например, число микрофонов может оказаться меньше числа точек сетки источника (в таком случае система (5) недоопределена); или шаг сетки источника во много раз меньше длины волны, что приводит к большим значениям числа обусловленности матрицы $\mathcal{T}_{a}^{{\text{*}}}{{\mathcal{T}}_{a}}$. Во многих работах по бимформингу в аналогичных случаях (см., например, [10]) предлагается использовать регуляризацию по Тихонову, то есть решать вместо исходной задачи задачу минимизации:
5.4.1. Недостаточное количество микрофонов. Рассмотрим регуляризацию на тестовой задаче из п.5.1 с количеством микрофонов M = 5, которое меньше количества точек сетки источника N = 10. Значения коэффициентов ${{a}_{S}} = 0.6$, ${{a}_{D}} = 0.12$, ${{a}_{{SD}}} = 5$ выбираем теми же самыми. Типичная зависимость C – нормы погрешности решения от параметра регуляризации α представлена на фиг. 9а. Аналогичный график для теста на поверхностях с параметрами $M = {{5}^{2}}$, $N = {{10}^{2}}$, ${{a}_{S}} = 0.8$, ${{a}_{D}} = 0.1$, ${{a}_{{SD}}} = 5$ приведен на фиг. 9б.
Фиг. 9.
Недостаточное число микрофонов. C-норма ошибки ${{\varepsilon }_{{1D}}}$ (а) и ${{\varepsilon }_{{2D}}}$ (б) в зависимости от параметра регуляризации α.
Видно, что в широком диапазоне значений параметра регуляризации $\alpha \in [{{10}^{{ - 4}}},{{10}^{{ - 2}}}]$ значениe ошибки не превышает нескольких процентов: ${{\varepsilon }_{{1D}}} = 0.02{\kern 1pt} - {\kern 1pt} 0.04$ и ${{\varepsilon }_{{2D}}} = 0.05{\kern 1pt} - {\kern 1pt} 0.06$.
5.4.2. Мелкая сетка источника. Необходимость регуляризации возникает также при использовании чересчур мелкой сетки источника с шагом distS существенно меньшими дифракционного предела. В качестве примера рассмотрим тестовую задачу из п. 5.1 с параметрами M = 15, N = 10, aSD = 5, aD = 0.12, но с тестовой функцией в виде гауссиана с полушириной 0.5λ и aS = 0.1. Зависимость нормы погрешности решения от параметра регуляризации α представлена на фиг. 10а. На фиг. 10б показана зависимость для аналогичного теста на поверхностях с M = 152, N = 102, aSD = 5, aS = 0.1, aD = 0.1. Отметим, что расчеты с α = 0, что означает отсутствие регуляризации, приводят к неприемлемо большой погрешности.
Фиг. 10.
Зависимость C-нормы ошибки ${{\varepsilon }_{{1D}}}$ (а) и ${{\varepsilon }_{{2D}}}$ (б) от параметра регуляризации α при значении ${\text{dis}}{{{\text{t}}}_{S}}$ много меньше дифракционного предела.
Из проведенных тестовых исследований можно сделать вывод, что значения ошибки не превышают ${{\varepsilon }_{{1D}}} = 0.025{\kern 1pt} - {\kern 1pt} 0.04$ в диапазоне $\alpha \in [{{10}^{{ - 3}}},{{10}^{{ - 2}}}]$ для одномерного случая (линия источника) и ${{\varepsilon }_{{2D}}} = 0.04{\kern 1pt} - {\kern 1pt} 0.08$ в диапазоне $\alpha \in [{{10}^{{ - 4}}},2 \times {{10}^{{ - 3}}}]$ для двумерного случая (поверхность источника). На фиг. 11 представлены карты амплитуд численного решения с регуляризацией при $\alpha = {{10}^{{ - 3}}}$ (а) и точного решения (б).
6. ТЕСТИРОВАНИЕ МЕТОДА БИМФОРМИНГА НА СИНТЕТИЧЕСКИХ ДАННЫХ
Численные тесты на корректность разработанного метода проведены по аналогии с [12], где для некоторой полосы $S$ носителя функции источника синтетические данные используются при тестировании решетки микрофонов в экспериментальном бимформинге. В отличие от этой работы, мы располагаем микрофоны на поверхности Фокса Вилльямса–Хокингса (ФВХ). Эти поверхности используются в вычислительном эксперименте по моделированию акустических полей, генерируемых сегментом крыла самолета с выпущенной механизацией на режиме посадки, описанном в разд. 7. Для наших тестов мы выбрали две такие поверхности D1 и D2 шириной 0.111 (см. фиг. 12). Полоса S для задания источника имеет длину 0.7 по оси X (вдоль хорды крыла) и ширину 0.06 по оси Z (по поперечному направлению).
Фиг. 12.
Геометрическая конфигурация тестовых задач с линией и поверхностью источника (рисунки (а) и (б) соответственно): линия и поверхность источника – S, ближняя и дальняя поверхности микрофонов – D1 и D2 соответственно.
6.1. Тесты с линией источника при низкой частоте
Зафиксируем длину звуковой волны λ = 0.1 и ${{\operatorname{M} }_{\infty }} = 0$. Поскольку λ > 0.06 – ширины поверхности S, мы не можем располагать более одной точки вдоль оси Z в соответствии с ограничением на шаг сетки функции источника ${\text{dis}}{{{\text{t}}}_{S}}\sim \lambda $. По оси X длина S достаточна, чтобы разместить несколько точек сетки источника. Отсюда и возникает постановка задачи бимформинга не для поверхности, а для линии источника. Зададим источник на отрезке $S = [0.15,{\text{\;}}0.85]$ кусочно-линейной функцией по точкам сетки ${{S}_{h}} = \{ 0.1,{\text{\;}}0.2,{\text{\;\;}}...,{\text{\;}}0.8\} $ на S со значениями $1$ в трех узлах X = 0.4, 0.5, 0.6 и 0 в остальных (см. фиг. 13а). Микрофоны располагаются на линии D1, их число M = 1663, и расположены они на поверхности ФВХ в узлах сетки, используемой в вычислительном эксперименте по моделированию турбулентного течения вокруг сегмента крыла с механизацией. Заметим, что для тестовых расчетов микрофоны могли бы располагаться на более простых поверхностях, но выбор поверхностей ФВХ сделан с целью единообразия с геометрией задачи из разд. 7. В табл. 1 для трех экспериментов на равномерных сетках ${{S}_{{h,1}}}$, ${{S}_{{h,2}}}$, ${{S}_{{h,3}}}$ функции источника с узлами на концах S указаны количество точек сетки N, шаги ${\text{dis}}{{{\text{t}}}_{S}}$, числа обусловленности ${{C}_{a}}~$ и L2-нормы погрешности. Соответствующие графики решений приведены на фиг. 13а.
Фиг. 13.
Точное решение и решения, полученные на разных сетках источника (а). Зависимость решения от числа узлов квадратуры на сетке ${{S}_{{h,1}}}$ (б).
Таблица 1.
Погрешность и число обусловленности для трех сеток источника
| Сетка источника | $N$ | ${\text{dis}}{{{\text{t}}}_{S}}$ | ${{C}_{a}}$ | ${{\left\| \varepsilon \right\|}_{{{{\mathcal{L}}_{2}}}}}$ |
|---|---|---|---|---|
| ${{S}_{{h,1}}}$ | 8 | 0.1 | 17 | 0.072 |
| ${{S}_{{h,2}}}$ | 15 | 0.05 | 22 | 0.00056 |
| ${{S}_{{h,3}}}$ | 29 | 0.025 | 6.5e+07 | 0.0031 |
Из графиков видно, что решения на сетках источника ${{S}_{{h,2}}}$ и ${{S}_{{h,3}}}$ отлично аппроксимируют точное решение. Заметное различие для грубой сетки источника ${{S}_{{h,1}}}$ (зеленая линия), которая имеет такой же шаг, что и сетка ${{S}_{h}}$ точной функции источника, связано только с тем, что их узлы сдвинуты на полшага друг от друга, в результате чего линейное восполнение решения на этих двух сетках дает разные функции. Заметим, что в дополнительном эксперименте с совпадением узлов обеих сеток графики точного и полученного решений совпадают.
Также можно отметить взрывной характер роста числа обусловленности ${{C}_{a}}~$ при уменьшении шага сетки источника от ${\text{dis}}{{{\text{t}}}_{S}} = 0.5\lambda $ до ${\text{dis}}{{{\text{t}}}_{S}} = 0.25\lambda $. При этом нарушается условие корректности матрицы бимформинга (см. п. 5.2). Тем не менее численное решение все еще достаточно хорошее, хотя можно наблюдать осцилляции вблизи X = 0.5.
На этой же задаче мы показываем, что точность интегрирования при вычислении элементов матрицы бимформинга в (7) существенно влияет на решение. На фиг. 13б приведены графики решений, полученных на сетке ${{S}_{{h,1}}}$ при разном количестве узлов квадратур Гаусса k. Заметим, что в случае k = 1 подход сводится к методу из работы [10] (Generalized Inverse Beamforming). Хорошо видно, что при малых k на грубой сетке ${{S}_{{h,1}}}$ решение оказывается неудовлетворительным. В табл. 2 приведены погрешность решения и число обусловленности ${{C}_{a}}$ матрицы бимформинга для различного числа узлов квадратур.
6.2. Тесты с плоскостью источника при высокой частоте
В этой серии тестов фиксируем длину волны λ = 0.01. Она уже меньше ширины полосы $S = [0.15,0.85] \times [ - 0.03,0.03]$ по направлению оси Z, поэтому имеет смысл рассматривать двумерную сетку источника. В качестве тестовой функции источника используем функцию той же трапециевидной формы по X, как и прежде, но сжатую в 10 раз. Она задана на равномерной сетке с постоянным шагом 0.01, принимает значение 1 в трех узлах $X = 0.49,{\text{\;}}0.50,{\text{\;}}0.51$ сетки и постоянна вдоль Z (см. фиг. 14а).
Для численного бимформинга берутся сетки ${{S}_{{h,1}}}$, ${{S}_{{h,2}}}$, ${{S}_{{h,3}}}$ источника, равномерные по обоим направлениям с шагом ${\text{dis}}{{{\text{t}}}_{S}} = 0.01,{\text{\;\;}}0.0085,{\text{\;\;}}0.005$ и количеством точек сетки источника N = 852, 1162, 3243 соответственно. Для сетки микрофонов используется поверхность D1 с количеством микрофонов M = 83 150. Полученные решения на сетках ${{S}_{{h,1}}}$, ${{S}_{{h,2}}}$ неотличимы на глаз от тестовой функции, которая изображена на фиг. 14а. Однако решение на сетке ${{S}_{{h,3}}}$, представленное на фиг. 14б, заметно отличается от точного тестового решения. Причиной этого является нарушение условия дифракционного предела (${\text{dis}}{{{\text{t}}}_{S}} > 0.6\lambda $). В табл. 3 для этих трех численных экспериментов указаны значения количества точек сетки источника N, шага сетки ${\text{dis}}{{{\text{t}}}_{S}}$, числа обусловленности ${{C}_{a}}$ и ${{L}_{2}}$-нормы погрешности. Несмотря на то что сетка ${{S}_{{h,1}}}$ совпадает с сеткой заданной тестовой функции источника по X, возникает небольшая погрешность решения. Это связано со спецификой обработки решения в текущей версии алгоритма на верхней и нижней границах полосы (Z = –0.03, 0.03).
Таблица 3.
Погрешность и число обусловленности Ca для каждой из сеток
| Сетка источника | $N$ | ${\text{dis}}{{{\text{t}}}_{S}}$ | ${{C}_{a}}$ | ${{\left\| \varepsilon \right\|}_{{{{\mathcal{L}}_{2}}}}}$ |
|---|---|---|---|---|
| ${{S}_{{h,1}}}$ | $71 \times 12$ | 0.01 | 81 | 0.00054 |
| ${{S}_{{h,2}}}$ | $83 \times 14$ | 0.0085 | 111 | 0.0036 |
| ${{S}_{{h,3}}}$ | $141 \times 23$ | 0.005 | 2.1e+17 | 0.014 |
Аналогично одномерному случаю исследуется влияние точности интегрирования при вычислении элементов матрицы в (7) на точность результата численного бимформинга. На фиг. 15 представлены решения на сетке ${{S}_{{h,1}}}$, полученные для двух значений k = 1, 3 числа узлов кубатур Гаусса в треугольнике (решение для k = 7 совпадает с тестовым, представленным на фиг. 14а). Более подробную информацию содержит табл. 4. Видно, что значения $k \leqslant 3$ не обеспечивают необходимую точность восстановления источника. Оптимальным значением оказывается k = 7.
6.3. Тесты с регуляризацией матрицы бимформинга
Для проверки работы алгоритма бимформинга с регуляризацией (см. п. 5.4) рассмотрим задачу восстановления тестовой функции источника – гауссиана с полушириной $4\lambda $, $\lambda = 0.05$, ${{\operatorname{M} }_{\infty }} = 0$, заданной на поверхности источника, представляющую собой полосу $[0.15,{\text{\;}}0.85] \times [ - 0.03,{\text{\;}}0.03]$. Поверхности микрофонов D1 или D2 имеют размер по ширине [–0.055, 0.055] (см. фиг. 16).
Сетка источника равномерная и состоит из 51 × 45 узлов. Соответственно, коэффициенты ${{({{a}_{S}})}_{X}} = 0.28$, ${{({{a}_{S}})}_{Y}} = 0.027$ заметно меньше предела корректности ${{a}_{S}} = 0.6$. Сетки микрофонов также равномерные на D1 и D2 и варьируются по числу узлов.
В табл. 5 представлены результаты численного исследования. Параметр регуляризации во всех расчетах берется равным $\alpha = {{10}^{{ - 4}}}$. Отметим, что без регуляризации число обусловленности ${{C}_{a}}$ превышает значение ${{10}^{{20}}}$, а задача с 871 микрофонами просто не имела бы смысла, поскольку M < N.
Таблица 5.
Параметры и погрешность решения трех расчетов с регуляризацией
| Сетка микрофонов | $M$ | ${{\left\| \varepsilon \right\|}_{{{{\mathcal{L}}_{2}}}}}$ |
|---|---|---|
| ${{D}_{2}}$ | 33 950 | 0.0072 |
| ${{D}_{1}}$ | 15 612 | 0.0069 |
| ${{D}_{1}}$ | 871 | 0.0070 |
Решения для первой и третьей тестовых задач представлены на фиг. 17. Отметим наличие пограничных эффектов для верхней и нижней границ полосы (как уже отмечалось, они возникают из-за специфики обработки решения на границах).
На фиг. 18 представлены графики точного и численных тестовых решений по линии Z = 0 полосы расположения источника. Можно сделать вывод, что регуляризация позволяет получать решение с высокой точностью даже при существенном нарушении условий построения корректной матрицы бимформинга (8), (10)–(13).
7. ОПРЕДЕЛЕНИЕ ЛОКАЛИЗАЦИИ И МОЩНОСТИ АКУСТИЧЕСКОГО ИСТОЧНИКА МЕТОДОМ БИМФОРМИНГА НА ОСНОВЕ АНАЛИЗА ДАННЫХ ВЫЧИСЛИТЕЛЬНОГО ЭКСПЕРИМЕНТА
Разработанный метод численного бимформинга был применен к пространственно-временным данным, полученным в ходе вычислительного эксперимента [16], который состоял в численном моделировании нестационарного турбулентного течения, генерируемого при обтекании аэродинамической модели 30P30N (прямое крыло самолета с выпущенной механизацией), и создаваемого при этом дальнего акустического поля. Расчет выполнялся для сегмента модели с поперечным размером 1/9 хорды крыла при угле атаки 5.5°. Характерное число Маха составляло ${{\operatorname{M} }_{\infty }} = 0.17$ при числе Рейнольдса Re = 1.7 × 106, посчитанному по длине хорды крыла. Моделирование турбулентного течения в [16] проводилось с использованием гибридного вихреразрешающего подхода IDDES (Improved Delayed Detached Eddy Simulation / Улучшенное моделирование отсоединенных вихрей с запаздыванием [17]) на неструктурированной гибридной сетке размером 35 миллионов узлов с использованием программного комплекса NOISEtte [18]. Акустическое излучение в дальнем поле, генерируемое вследствие взаимодействия набегающего турбулентного течения и крылового сегмента с механизацией, рассчитывалось на стадии постпроцессинга данных согласно интегральному методу ФВХ.
Данные экспериментального исследования крыловой конфигурации 30P30N представлены в работе [19], а результаты численных исследований, выполненных различными научными коллективами с помощью известных программных продуктов и исследовательских кодов, – в статье [20].
На фиг. 19 изображено мгновенное распределение поля давления, формирующегося при нестационарном дозвуковом обтекании сегмента трехкомпонентного крыла с механизацией. Можно заметить наличие источников звуковых волн разной длины в окрестности зазора между предкрылком и крылом, зазора между крылом и закрылком, а также на задней кромке закрылка.
Фиг. 19.
Полученное в вычислительном эксперименте мгновенное поле давления, формирующееся при турбулентном обтекании крыла с механизацией.
В качестве примера будем исследовать распределенный акустический источник, формирующийся в зазоре между предкрылком и крылом. Спектральная плотность мощности пульсаций давления в выбранной точке Р8 (фиг. 20а) в окрестности предкрылка приведена на фиг. 20б. В представленном спектре с выделим две третьоктавные полосы частот с центрами при числах Струхаля Sh = 21 и Sh = 180 (безразмерное число Струхаля определяется как ${\text{Sh}} = \frac{{fL}}{V}$, где f – частота, V – скорость, L – характерный размер). Выберем среди различимых пиков в спектре (фиг. 20б) пик, соответствующий числу Струхаля Sh = 21. Второй пик, интересный для теста разработанного подхода, соответствует Sh = 180 и находится в высокочастотной части спектра за пределами графика.
Фиг. 20.
Мгновенное поле давления вблизи предкрылка (а) и спектральная плотность мощности акустического излучения в точке P8 (б).
Полученные в ходе вычислительного эксперимента [16] данные накапливались на окружающих всю конфигурацию контрольных ФВХ поверхностях. Мы используем две из них в качестве поверхностей D1 или D2, где расположены виртуальные микрофоны (фиг. 21–23а).
Фиг. 21.
Трехзвенная линия источника S и две поверхности микрофонов (а). Распределение амплитуды низкочастотного источника (б): на линии (1D) и вдоль средней линии на кусочно-плоской поверхности (2D).
Фиг. 22.
Кусочно-плоская поверхность низкочастотной функции источника и поверхность микрофонов (а). Распределение амплитуды низкочастотного источника на поверхности при использовании регуляризации (б).
Фиг. 23.
Кусочно-плоская поверхность высокочастотного источника и поверхность микрофонов (а). Распределение амплитуды высокочастотного источника на поверхности (б).
Полученные временные ряды, представляющие собой пульсации давления в точках сетки виртуальных микрофонов, подвергались следующей обработке:
1) применялись высокочастотный и низкочастотный фильтры, вырезающие интересующие нас полосы частот;
2) сигналы микрофонов получались путем усреднения по случайным временным окнам из всего отрезка времени счета.
Для проведения усреднения выбиралось достаточно большое количество так называемых “сценариев”, каждый из которых характеризовался случайно заданным начальным моментом времени и несколькими последовательными временными окнами заданной длины. Для каждого из окон вычислялся спектр и проводилось усреднение всех спектров по окнам в сценарии. Затем полученные таким образом спектры усреднялись по всем сценариям. В результате были получены устойчивые к изменениям параметров сценариев спектральные характеристики.
7.1. Источник низких частот в окрестности числа Струхаля 21
Рассматривается интервал частот $18.7 < \operatorname{Sh} < 23.6$ с центральной частотой ${\text{S}}{{{\text{h}}}_{c}} = 21$ и длиной волны $\lambda = c{\text{Sh}}_{c}^{{ - 1}} = \operatorname{M} _{\infty }^{{ - 1}}{\text{Sh}}_{c}^{{ - 1}} \approx 0.28$. Из-за небольшой ширины сегмента крыла $H \approx 0.111$, соответствующей 1/9 длины хорды, характерная длина волн в рассматриваемом диапазоне частот $\lambda > H$, а потому функция источника может представлять собой только константу вдоль оси Z. Поэтому будем рассматривать задачу для линии расположения источника, заметив, что при этом данные для микрофонов предварительно усредняются по ширине крылового сегмента.
Линию источника задаем в виде ломаной (фиг. 21) и строим равномерную сетку на каждом из трех ее звеньев с шагом ${\text{dis}}{{{\text{t}}}_{S}} \approx 0.24\lambda $. Это значение находится ниже дифракционного предела ($0.4\lambda $), но все еще может использоваться, что, в частности, подтверждается сравнением с результатом расчета с регуляризацией, который производился для случая расположения источника на поверхности (см. фиг. 22) с двумерной сеткой источника размером 21× 5 при параметре регуляризации $\alpha = {{10}^{{ - 4}}}$.
Выполненные расчеты и, в частности, анализ распределения амплитуды (фиг. 21б, фиг. 22) говорят о существовании достаточно мощного акустического источника, излучающего в выбранном низкочастотном диапазоне волн в окрестности Sh = 21, при X ≈ 0.1, что соответствует зазору между предкрылком и крылом.
7.2. Источник высоких частот в окрестности числа Струхаля 180
Рассматривается интервал чисел Струхаля $160 < \operatorname{Sh} < 202$ с центральной частотой ${\text{S}}{{{\text{h}}}_{c}} = 180$ и длиной волны $\lambda \approx 0.033$. Акустический источник расположен на поверхности, состоящей из трех плоских участков, на сетке размером 61 × 5. Согласно полученным в работе критериям корректности, матрица бимформинга при данных параметрах корректная. Результаты нахождения источника в заданной полосе высоких частот показаны на фиг. 23.
Выполненные расчеты и, в частности, анализ распределения амплитуды (фиг. 23а) выявляют акустический источник, излучающий в выбранном высокочастотном диапазоне волн в окрестности Sh = 180, при X ≈ 0.1, что также соответствует зазору между предкрылком и крылом. Вместе с тем можно видеть, что высокочастотный источник по мощности существенно слабее ранее выявленного источника низкой частоты. Более того, пульсации давления в столь высоком диапазоне частот вполне могут быть обусловлены нефизическим шумом численного происхождения. Однако эти слабые пульсации присутствовали в данных вычислительного эксперимента и разработанный метод бимформинга вполне четко их улавливает.
8. ЗАКЛЮЧЕНИЕ
В работе разработан метод численного бимформинга для определения непрерывной (кусочно-линейной) функции акустического источника на основе анализа пространственно-временных данных вычислительного эксперимента по моделированию нестационарного турбулентного течения. Естественным условием успешности такого подхода является высокая точность входных данных, которыми являются воспроизводимые в расчете параметры течения и его пульсационные характеристики. Поэтому предполагается, что моделирование турбулентного течения проводится с использованием вихреразрешающих подходов, численных алгоритмов повышенной точности, а также пространственных сеток большой размерности. Потребность в значительных вычислительных ресурсах для таких расчетов обеспечивается, в свою очередь, современными суперкомпьютерными системами с растущей производительностью.
В статье выписана математическая постановка задачи бимформинга с переходом из временной в частотную область, рассмотрен подход к ее дискретизации и проведена верификация разработанного численного метода на тестовых экспериментах с синтетическими данными. Особенностью предложенной дискретизации является использование непрерывных базисных функций при аппроксимации распределенного источника в сочетании с кубатурными формулами высокого порядка для вычисления элементов матрицы бимформинга.
Исходя из физических соображений и требований аппроксимации, сформулированы условия на параметры сетки дискретизации источника монопольного типа и сетки микрофонов для обеспечения корректности и точности работы алгоритма. Предложен способ оценки входящих в указанные условия численных параметров, а именно, шагов сетки функции источника и микрофонов, расстояния между ними, а также число узлов в кубатурных формулах Гаусса. Корректность задания этих параметров обеспечивается путем контроля числа обусловленности оператора бимформинга. Так, резкий рост числа обусловленности наблюдается, если использовать: а) шаг сетки аппроксимации источника ниже некоторой величины, соизмеримой с длиной волны рассматриваемой частоты; б) недостаточное число узлов в кубатурных формулах Гаусса. Показано, что при локализации источника монопольного типа на линии необходимо, чтобы отношение шага сетки источника к длине волны было не менее 0.4, а число узлов в кубатурных формулах Гаусса – не менее 4. Аналогичные параметры для случая локализации источника на поверхности равны 0.6 и 7 соответственно. При нарушении условий корректности рассмотрен также подход регуляризации оператора бимформинга с целью сохранения возможности получения функции источника с приемлемой точностью.
Продемонстрировано, что предложенный метод численного бимформинга обеспечивает высокую точность восстановления функции источника на тестовых примерах, а его использование к данным вычислительного эксперимента по турбулентному обтеканию сегмента крыла с механизацией дает хорошее согласие с ранее полученными результатами анализа акустического поля.
Авторы благодарят своих коллег А.П. Дубеня, А.В. Горобца и П.В. Родионова за полезные дискуссии и замечания.
Список литературы
Johnson H., Dudgeon D.E. Array Signal Processing: Concepts and Techniques. London: Prentice-Hall, 1993. 533 p.
Dougherty R.P. Phased-array Beamforming Applied to Aeroacoustics // 15th AIAA/CEAS Aeroacoustics Conference (30th AIAA Aeroacoustics Conference), 2009.
Capon J. High resolution frequency-wavenumber spectrum analysis // Proc. of the IEEE. 1969. V. 57. № 8. P. 1408–1418.
Cox H., Zeskind R.M., Owen M.M. Robust adaptive beamforming // IEEE Trans. Acoust. Speech Signal Processing. 1987. V. 35. P. 1365–1376.
Högbom J.A. Aperture synthesis with a non-regular distribution of interferometer baselines // Astronomy and Astrophysics Supplement Series. 1974. V. 15. P. 417–426.
Brooks T.F., Humphreys W.M. A deconvolution approach for the mapping of acoustic sources (DAMAS) determined from phased microphone arrays // J. Sound Vib. 2006. V. 294. P. 856–879.
Dougherty R.P. Extension of DAMAS and benefits and limitations of deconvolution in beamforming // AIAA Paper 2005–2961. 2005.
Brooks T.F., Humphreys W.M. Extension of DAMAS phased array processing for spatial coherence determination(DAMAS-C) // AIAA-2006-2654. 2006.
Sijtsma P. CLEAN based on spatial source coherence // AIAA-2007-3436. 2007.
Suzuki T. Generalized Inverse Beam-forming Algorithm Resolving Coherent/Incoherent, Distributed and Multipole Sources // J. Sound Vib. 2011. V. 330. P. 5835–5851.
Merino-Martínez R., Sijtsma P., Snellen M., Ahlefeldt T., Antoni J., Bahr C.J., Blacodon D., Ernst D., Finez A., Funke S., Geyer T.F., Haxter S., Herold G., Huang X., Humphreys W.M., Leclère Q., Malgoezar A., Michel U., Padois T., Pereira A., Picard C., Sarradj E., Siller H., Simons D.G., Spehr C. A review of acoustic imaging methods using phased microphone arrays // CEAS Aeronautical J. 2019. V. 10. № 1. P. 197–230.
Sarradj E., Herold G., Sijtsma P., Merino-Martínez R., Geyer T.F., Bahr C.J., Porteous R., Moreau D., Doolan C.J. A Microphone Array Method Benchmarking Exercise using Synthesized Input Data // AIAA paper 2017–3719. 2017.
Голованов Н.Н. Геометрическое моделирование. М.: Физматлит, 2002. 472 с.
Ландау Л.Д., Лифшиц Е.М. Гидродинамика. Теоретическая физика: т. VI. 3-е изд., перераб. М.: Физматлит, 1986. 736 с.
Wu T.W., Lee L. A Direct Boundary Integral Formulation For Acoustic Radiation In A Subsonic Uniform Flow // J. Sound Vib. 1994. V. 175. P. 51–63.
Bobkov V., Gorobets A., Duben A., Kozubskaya T., Tsvetkova V. Towards affordable CAA simulations of airliner’s wings with deployed high-lift devices // Book of Abstracts of the 5-th Internat. Workshop “Computational Experiment in AeroAcoustics“. 2018. P. 36–37.
Shur M.L., Spalart P.R., Strelets M.Kh.,Travin A.K. A hybrid RANS-LES approach with delayed-DES and wall-modeled LES capabilities // Internat. Journal of Heat and Fluid Flow. 2008. V. 29. № 6. P. 1638–1649.
Абалакин И.В., Бахвалов П.А., Горобец А.В., Дубень А.П., Козубская Т.К. Параллельный программный комплекс NOISETTE для крупномасштабных расчетов задач аэродинамики и аэроакустики // Вычисл. методы и программирование. 2012. Т. 13. С. 110–125.
Pascioni K., Cattafesta L.N., Choudhari M.M. An Experimental Investigation of the 30P30N Multi-Element High-Lift Airfoil // AIAA paper 2014–3062. 2014.
Choudhari M.M., Lockard D.P. Assessment of Slat Noise Predictions for 30P30N High-Lift Configuration from BANC-III Workshop // AIAA paper 2015–2844. 2015.
Дополнительные материалы отсутствуют.
Инструменты
Журнал вычислительной математики и математической физики
Пн.-пт.: 10.00-19.00