Application of spectral proper orthogonal decomposition to the analysis of sound field of aeroacoustic sources

Abstract

Spectral proper orthogonal decomposition (SPOD) is proposed for the identification of the multipole structure of aeroacoustic noise sources from far-field measurements. The method is verified via tests with point multipoles and validated using experimental data on flow-induced cylinder noise and turbulent jet noise.

Full Text

Введение

В задачах аэроакустики часто приходится иметь дело с источниками различной мультипольности. Согласно одной из общих формулировок акустической аналогии [1], пульсации массового расхода соответствуют монопольному излучению, пульсации силы – дипольному, пульсации рейнольдсовских напряжений – квадрупольному. В некоторых течениях значимыми могут быть и мультиполи более высокого порядка [2]. Таким образом, идентификация различных мультипольных компонент шума, излучаемого исследуемым источником, может быть полезна с точки зрения бесконтактной диагностики физической природы самого источника. Такая идентификация не всегда возможна, поскольку в случае произвольного пространственного распределения источников различной мультипольности задача определения их типа, положения и интенсивности не имеет однозначного решения [3-5]. Тем не менее, в некоторых практически реализуемых случаях постановка подобной задачи может быть не лишена смысла. Так, для относительно медленных дозвуковых течений аэроакустические источники часто являются компактными [3] (во всех или в выделенных направлениях), что существенно упрощает задачу их исследования. Примерами компактных источников являются вихревые кольца [6, 7], тела, обтекаемые дозвуковым потоком [8-10], турбулентные дозвуковые струи [2, 11, 12], плазменные актуаторы на умеренных частотах [13, 14].

При этом, проводя анализ мультипольной структуры акустического излучения сложных источников, следует различать эффекты, связанные с процессами генерации шума, и эффекты, обусловленные особенностями распространения шума. Как известно [15], решение волнового уравнения в безграничной неподвижной среде на сферической поверхности, окружающей источник звука, может быть представлено в виде разложения по полной системе сферических функций $Y_l^n=e^{in\phi }{P}_l^{\left|n\right|}(\cos \theta )$, где $\phi$ и θ – соответственно азимутальный и полярный углы сферической системы координат $(r,\theta ,\phi )$, $P_l^{\left|n\right|}(\cos \theta )$ – присоединенные полиномы Лежандра, $l\ge \left|n\right|$. Для акустически компактных неподвижных источников, расположенных в области начала сферической системы координат, имеет место однозначное соответствие между его “физическими” мультипольными компонентами, связанными с теми или иными механизмами генерации шума, и сферическими функциями, т.е. значение l соответствует порядку физического мультиполя: l = 0 – монопольное излучение, содержащее только осесимметричную азимутальную компоненту n = 0; l = 1 – дипольное, состоящее из двух азимутальных компонент |n| = 0, 1; l = 2 – квадрупольное, содержащее три азимутальные компоненты |n| =0, 1, 2 и т.д. [9, 16].

Такая однозначная связь нарушается в случае, если источник смещен относительно начала координат [9, 17, 18] (при этом смещение, разумеется, должно быть сравнимо или больше длины излучаемой волны), повернут [19], находится в движении [2] (эффекты конвекции) и/или излучает в неоднородную движущуюся среду [20, 21] (эффекты рефракции). В этом случае представление направленностей физических мультиполей через исходные сферические функции, хотя и остается справедливым, перестает быть достаточно наглядным в тем большей степени, чем сильнее выражены эффекты смещения, движения и поворота, поскольку физические мультиполи “распределяются” по различным сферическим гармоникам.

Часто в задачах аэроакустики указанные трансляционные эффекты (не включающие поворот) носят одномерный характер, т.е. имеется выделенное направление, вдоль которого источник может быть смещен, может двигаться и/или вдоль которого движется среда. Обычно это направление совпадает с направлением набегающего потока. В таком случае, если совместить полярную ось системы координат с этим выделенным направлением так, чтобы она при этом проходила через область источников, азимутальная структура звукового поля физических мультиполей не поменяется, а изменится лишь их полярная направленность. Отметим, что влияние на диаграмму направленности источника его сдвига вдоль выделенного направления всегда может быть сделано пренебрежимо малым при достаточном удалении от источника, поэтому главными эффектами распространения, которые могут исказить (по сравнению с излучением в покоящейся среде) вид физических мультиполей, являются эффекты конвекции и рефракции. Таким образом, наглядность мультипольного разложения с точки зрения его связи с физическими механизмами генерации шума можно сохранить, перейдя от системы функций $e^{in\phi }{P}_l^{\left|n\right|}(\cos \theta )$ к системе $e^{in\phi }{f}_l^{\left|n\right|}(\theta ,{\text{M}}_s,\text{M})$, где функция $f_l^{\left|n\right|}(\theta ,{\text{M}}_s,\text{M})$ описывает полярную диаграмму направленности базисных физических мультиполей с учетом эффектов конвекции/рефракции, величины ${\text{M}}_s$ и M обозначают числа Маха движения источника и среды соответственно, и выполняется $f_l^{\left|n\right|}(\theta ,{\text{M}}_s,\text{M})\to P_l^{\left|n\right|}(\cos \theta )$ при ${\text{M}}_s,\text{M}\to 0$. Это означает, что направленности отдельных азимутальных компонент звукового поля, определенных относительно оси, направленной вдоль среднего течения и проходящей через область источников, могут отражать особенности механизмов генерации шума, если их интерпретировать в указанном выше смысле.

Именно такая интерпретация азимутальных компонент звукового поля аэроакустических источников легла в основу метода азимутальной декомпозиции (МАД), разработанного в ЦАГИ [11, 22]. Вид “деформированных” эффектами распространения направленностей базисных мультиполей в некоторых случаях может быть оценен теоретически [2, 20, 21, 23, 24], что позволяет проводить детальное сравнение моделей источников шума и результатов измерений [2, 11, 20, 22, 24]. Таким образом, определение азимутального состава звукового поля аэроакустических источников является мощным инструментом анализа физических механизмов генерации шума.

Измерения азимутального состава звукового поля обычно проводятся с помощью азимутальных решеток микрофонов [11, 22, 25-27], причем число микрофонов в каждом сечении должно быть не меньше, чем число значимых азимутальных компонент, хотя в случае звуковых полей, излучаемых осесимметричными в среднем течениями, необходимое число микрофонов может быть уменьшено вдвое [28]. Альтернативным методом является использование метода азимутального разложения корреляционной функции, для чего достаточно пары микрофонов, один из которых неподвижен и является опорным, а второй – может двигаться по азимутальной координате вокруг источника [29, 30]. В этом случае можно разрешить большее число азимутальных мод, чем при использовании нескольких неподвижных микрофонов, однако измерительная система оказывается более сложной для реализации из-за наличия подвижного датчика, а полученные результаты позволяют говорить лишь об энергетическом составе суммарного поля (фазовые характеристики мод не разрешаются).

В некоторых случаях, из-за ограничений экспериментальной установки и/или измерительной системы, ни один из указанных выше способов азимутальной декомпозиции не может быть реализован. Такие ограничения могут возникать и в лабораторных условиях из-за различных геометрических ограничений (загромождение оборудованием, близость отражающих поверхностей и др.), и, например, при акустических измерениях натурных объектов на больших установках [28], когда шум в дальнем поле можно измерять только под заданным азимутальным углом или в узком азимутальном секторе. В таком случае становится актуальным вопрос о том, возможно ли оценить мультипольную структуру источника, используя измерительную систему, расположенную в одной азимутальной плоскости (по сути – обычный массив микрофонов, распределенных по дуге окружности вокруг источника).

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

В настоящей работе для такой оценки предлагается использовать метод собственного ортогонального разложения в частотной области (так называемый метод SPOD – spectral proper orthogonal decomposition). Метод SPOD в последнее десятилетие занял прочное место в ряду инструментов, использующихся для анализа турбулентных течений, позволяя, например, идентифицировать в таких течениях пространственно-временные когерентные структуры, проводить детальное сопоставление теоретических моделей с результатами физических или численных экспериментов, строить модели пониженного порядка для сложных течений с целью реализации активного управления ими (см., например, обзор [34] и соответствующие ссылки в нем). Детальное описание метода SPOD приведено в [35], где также обсуждается его связь со ставшим уже классическим в гидродинамике методом собственного ортогонального разложения (POD – proper orthogonal decomposition) [36], который, в свою очередь, восходит к методу главных компонент, предложенному Пирсоном [37]. Применительно к задачам аэроакустики метод SPOD зачастую используется для выделения когерентных структур в турбулентной области течения и, в некоторых случаях, последующей оценке звукового излучения от таких структур [38, 39]. В настоящей работе, в отличие от указанных исследований, предлагается построить метод SPOD для анализа структуры излучения акустически компактных источников в дальнем поле так, чтобы в результате измерений можно было сделать вывод о типе мультипольного источника и оценить интенсивности его отдельных мультипольных компонент.

Работа структурирована следующим образом. В разделе 1 описана идея метода на базе SPOD-разложения и оценка его работоспособности на модельных примерах, в разделе 2 представлено описание измерений, проведенных в заглушенной камере АК-2 ЦАГИ и результатов применения разработанного метода к анализу полученных экспериментальных данных. В Заключении сформулированы основные результаты работы.

1. Оценка возможности применения spod к анализу структуры звукового поля компактных источников

Приведем краткое описание метода SPOD, следуя [35]. Пусть $u(x,t)$ – одномерное стационарное случайное поле с нулевым средним. Во многих задачах, особенно в тех, где имеет место пространственный перенос возмущений, удобно оперировать со спектральным представлением данного поля в виде $\widehat{u}(x,t)$ и следить за поведением отдельных гармоник, сравнивая его с различными теоретическими моделями.

Метод SPOD позволяет найти оптимальное разложение случайной величины $\widehat{u}(x,f)$ по системе детерминированных ортогональных функций ${\psi }_j(x,f)$, упорядоченных по убыванию энергии, так что

 $\widehat{u}(x,f)={\displaystyle \sum _{j=1}^{\infty }{a}_j(f){\psi }_j(x,f)}$,       (1)

причем оптимальность понимается в том смысле, что при взятии N первых членов разложения (1) точность приближения $\widehat{u}(x,f)$ для любого другого набора базисных функций будет ниже. Для формулировки и решения задачи вводится скалярное произведение

${⟨u,\nu ⟩}_{x,t}={\displaystyle \iint {\nu }^{\ast }(x,t)W(x)u(x,t)dxdt}$,    (2)

определяющее понятие нормы и ортогональности рассматриваемых функций. Здесь * обозначает комплексное сопряжение (и транспонирование – для матриц), $W(x)$ – неотрицательная весовая функция, выбирающаяся из физических соображений, как, например, это будет показано ниже. Для стационарных случайных полей все утверждения об оптимальности, наилучшем приближении и т.д. понимаются в среднем, т.е. в смысле математического ожидания соответствующих случайных величин.

Можно показать [35, 36], что в частотной области задача определения функций ${\psi }_j(x,f)$ сводится к следующей задаче на собственные значения

 ${\displaystyle \int S(x,x^{\prime },f)W(x^{\prime })\psi (x^{\prime },f)d{x}^{\prime }}=\lambda (f)\psi (x,f)$,   (3)

где $S(x,x^{\prime },f)={\displaystyle \int C(x,x^{\prime },\tau )e^{-i2\pi f\tau }d\tau }$ – Фурье-образ тензора взаимных корреляций $C(x,x^{\prime },t-t^{\prime })=E\{u(x,t)u^{\ast }(x^{\prime },t^{\prime })\}$ $C(x,x^{\prime },t-t^{\prime })=E\{u(x,t)u^{\ast }(x^{\prime },t^{\prime })\}$, где $E\left\{g\right\}$ обозначает математическое ожидание случайной величины g и использована гипотеза о стационарности рассматриваемых случайных процессов. Можно показать, что в указанной постановке для задачи (3) существует счетное множество собственных мод, т.е. пар $\{{\psi }_j,{\lambda }_j\}$, которые могут быть ранжированы по величине (неотрицательных) собственных чисел ${\lambda }_1\ge {\lambda }_2\ge \mathrm{...}\ge 0$. Собственные функции ${\psi }_j(x,f)$ – SPOD-моды – являются взаимно ортогональными в смысле пространственного скалярного произведения

 ${⟨{\psi }_j,{\psi }_k⟩}_x={\displaystyle \int {\psi }_k^{\ast }(x,f)W(x){\psi }_j(x,f)dx}={\delta }_{jk}$, (4)

где ${\delta }_{jk}$ – дельта Кронекера. Кроме того, все базисные гармоники ${\psi }_j(x,f)e^{i2\pi ft}$ взаимно ортогональны и в смысле исходного пространственно-временного скалярного произведения (2). Таким образом, спектр $\widehat{u}(x,f)$ каждой отдельной реализации рассматриваемого случайного поля может быть оптимальным образом представлен в виде (1), где коэффициенты $a_j={⟨\widehat{u},{\psi }_j⟩}_x$ являются проекциями исходного поля на базисные функции, причем эти коэффициенты взаимно некоррелированы и $E\left\{a_j{a}_k^{\ast }\right\}={\lambda }_j{\delta }_{jk}$, а значит, взаимная спектральная плотность $S(x,x^{\prime },f)$ представляется в диагональной форме

 $S(x,x^{\prime },f)={\displaystyle \sum _{j=1}^{\infty }{\lambda }_j(f){\psi }_j(x,f){\psi }_j^{\ast }(x^{\prime },f)}$.             (5)

Такого рода представление случайного поля будет полезно в том случае, если в выражение (5) вносит вклад малое число мод. Это позволяет строить низкопорядковые модели исходного явления, а также пытаться сопоставить отдельным модам конкретные физические механизмы.

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

Пусть в начале координат расположены некоррелированные мультиполи порядка l. Тогда дальнее звуковое поле может быть представлено в виде

$\begin{array}{c}p(r,\theta ,\phi ,t)\approx \frac{1}{r}{\displaystyle \sum _{n=0}^l\left(a_n(t-r/c)\cos n\phi +\right.}\\ \left.+b_n(t-r/c)\sin n\phi \right)P_l^n(\cos \theta ),\end{array}$               (6)

где функции $a_n$ и $b_n$ взаимно некоррелированы. В частотной области (6) приобретает вид

$\widehat{p}(r,\theta ,\phi ,f)\approx \frac{e^{-i2\pi fr/c}}{r}{\displaystyle \sum _{n=0}^l({\widehat{a}}_n(f)\cos n\phi +{\widehat{b}}_n(f)\sin n\phi P_l^n(\cos \theta )}.$ Для точек наблюдения, расположенных в одной азимутальной плоскости $\phi =const$, без потери общности можно написать

 $\begin{array}{c}\widehat{p}(r,\theta \mathrm{,0,}f)\equiv \tilde{p}(r,\theta ,f)\approx \\ \approx \frac{e^{-i2\pi fr/c}}{r}{\displaystyle \sum _{n=0}^l{\widehat{a}}_n(f)P_l^n(\cos \theta )}.\end{array}$    (7)

Взаимная спектральная плотность для двух любых точек в этой плоскости в дальнем поле, с учетом взаимной некоррелированности отдельных мультиполей, будет выражаться как

$\begin{array}{c}\widehat{S}(r,\theta ,r^{\prime },{\theta }^{\prime },f)\approx \\ \approx {\displaystyle \sum _{n=0}^l{A}_n(f)P_l^n(\cos \theta )\frac{e^{-i2\pi fr/c}}{r}{P}_l^n(\cos {\theta }^{\prime })\frac{e^{i2\pi f{r}^{\prime }/c}}{r^{\prime }},}\end{array}$            (8)

где $A_n$ соответствует интенсивности источников, т.е. $E\left\{{\widehat{a}}_n{\widehat{a}}_m^{\ast }\right\}=A_n{\delta }_{nm}$. Отметим сходство выражений (8) и (5). Рассмотрим, каким образом можно поставить задачу SPOD-разложения звукового поля в заданной азимутальной плоскости, чтобы SPOD-моды были близки функциям $P_l^n(\cos \theta )e^{-i2\pi fr/c}/r$, характеризующим источник шума. Для простоты положим, что точки наблюдения лежат на дуге окружности $r=R$, как это обычно имеет место в акустических экспериментах. В этом случае (8) сводится к

$S(\theta ,{\theta }^{\prime },f)\approx {\displaystyle \sum _{n=0}^l{A}_n(f)P_l^n(\cos \theta )P_l^n(\cos {\theta }^{\prime })}$.           (9)

В соответствии с (4), (5) функции $P_l^n(\cos \theta )$ будут SPOD-модами, если для них будет выполняться условие ортогональности

 ${⟨P_l^n,P_l^m⟩}_{\theta }={\displaystyle \underset{0}{\overset{\pi }{\int }}{P}_l^m(\cos \theta )W(\theta )P_l^n(\cos \theta )d\theta }={\delta }_{nm}$. (10)

Как известно, присоединенные полиномы Лежандра не являются взаимно ортогональными, ортогональны только некоторые их подмножества [40]. В рассматриваемой задаче мы интересуемся возможностью выделения отдельных компонент мультиполей заданной степени l. Для такого подмножества $l=const$ имеет место ортогональность вида

${\displaystyle \underset{-1}{\overset{1}{\int }}{P}_l^m(q)P_l^n(q)\frac{dq}{1-q^2}}=\left[\begin{array}{l}{\delta }_{mn}\frac{(l+n)!}{(l-n)!n},n>\mathrm{0,}\\ \infty ,n=\mathrm{0,}\end{array}\right.$     (11)

где $q=\cos \theta$ и формально показано, что при $n=0$ соответствующий интеграл расходится. В строгой постановке расходящийся интеграл может быть проблемой, однако мы далее сосредоточимся на случаях, реализуемых в физических или численных экспериментах, в которых приходится иметь дело с дискретным набором координатных точек, в нашем случае – углов наблюдения $\left\{{\theta }_i\right\}$, $i=\mathrm{1,...,}N$. Положим, что углы распределены равномерно с шагом $\Delta \theta$ в секторе $0<{\theta }_1<\mathrm{...}<{\theta }_N<\pi$. Тогда условие (11) можно записать в виде

 ${\displaystyle \sum _{i=1}^N{P}_l^m(\cos {\theta }_i)P_l^n(\cos {\theta }_i)\frac{\Delta \theta }{\sin {\theta }_i}}\approx {\delta }_{mn}{g}_{ln}$,                   (12)

где $g_{ln}>0$ – набор положительных чисел, которые при достаточно мелком разбиении отрезка $[\mathrm{0,}\pi ]$ близки к $(l+n)!/((l-n)!n)$ при $n>0$ и к ${\displaystyle \underset{\cos {\theta }_N}{\overset{\cos {\theta }_1}{\int }}{\left(P_l^0(q)\right)}^2{\left(1-q^2\right)}^{-1}dq}$ – при $n=0$. Если понимать условие ортогональности в смысле (12), т.е. определить весовую функцию как

$W(\theta )=\frac{\Delta \theta }{\sin \theta }$,                     (13)

то соответствующие SPOD-моды, ортогональные в смысле так сконструированного скалярного произведения

 ${⟨{\psi }_j,{\psi }_k⟩}_{\theta }={\displaystyle \sum _{i=1}^N{\psi }_k^{\ast }({\theta }_i,f)W({\theta }_i){\psi }_j({\theta }_i,f)}={\delta }_{jk}$, (14)

будут близки к соответствующим присоединенным полиномам Лежандра (с точностью до множителей $g_{ln}$ и ошибки, связанной с дискретизацией по углу $\theta$ ), а спектральная плотность излучения в каждом направлении

$S(\theta ,\theta ,f)={\displaystyle \sum _{j=1}^N{\lambda }_j(f)|{\psi }_j(\theta ,f){|}^2}$                               (15)

будет, таким образом, разложена по базисным мультипольным компонентам в порядке убывания энергии. Отметим, что при этом неважно, будут ли коэффициенты $g_{ln}$ отнесены к собственному числу или к собственной функции, поскольку интенсивность излучения каждого мультиполя характеризуется произведением ${\lambda }_j(f)|{\psi }_j(\theta ,f){|}^2$. Число значимых собственных чисел в таком SPOD-разложении должно быть равно l+1 – числу базисных мультиполей заданной степени, дающих вклад в излучение в фиксированной азимутальной плоскости.

Задача (3), решением которой являются собственные числа и собственные функции, в дискретной форме для каждой частотной компоненты f принимает следующий вид [35]

 $S_{f}W{\Psi }_f={\Psi }_f{\Lambda }_f$,           (16)

где все матрицы имеют размер $N\times N$ и состоят из следующих элементов: ${\Psi }_{fij}={\psi }_j({\theta }_i,f)$, ${\Lambda }_{fij}={\lambda }_j(f){\delta }_{ij}$, $W_{ij}=\Delta \theta /\sin {\theta }_i{\delta }_{ij}$, $S_{fij}$ – элементы матрицы взаимных спектров для частоты f. Для определения $S_{fij}$ используется стандартный метод Уэлча [41] в виде, предложенном в [35] и кратко описанном ниже.

Пусть общая длина реализации случайного процесса, зафиксированного в точках наблюдения, составляет $M$ дискретных отсчетов с шагом по времени $\Delta t$. Соответствующий набор сигналов для $N$ точек наблюдения удобно сформировать в матрицу $Q=[q_1\mathrm{,...,}{q}_M]$ размером $N\times M$, в строках которой записаны временные реализации в соответствующих точках наблюдения, а в столбцах – мгновенные “слепки” случайного поля на наборе точек наблюдения. Далее полная временная реализация разбивается на $N_b$ блоков длиной $N_f$ каждый, и для каждого блока с номером n формируется матрица $Q^{(n)}=[q_1^{(n)}\mathrm{,...,}{q}_{N_f}^{(n)}]$ размером $N\times N_f$, состоящая из соответствующих столбцов матрицы $Q$. Затем матрицы $Q^{(n)}$ переводятся в частотное пространство путем применения к их строкам дискретного преобразования Фурье ${\widehat{q}}_k^{(n)}=N_f^{-1/2}{\displaystyle {\sum }_{j=1}^{N_f}{q}_j^{(n)}\exp (-i2\pi f_k(j-1)\Delta t)}$, где $f_k=(k-1)/(N_f\Delta t)$, $k=\mathrm{1,...,}{N}_f$ (при необходимости, преобразование может выполняться с весовой функцией), в результате чего формируются соответствующие $Q^{(n)}$ матрицы ${\widehat{Q}}^{(n)}=[{\widehat{q}}_1^{(n)}\mathrm{,...,}{\widehat{q}}_{N_f}^{(n)}]$. Тогда взаимная спектральная матрица оценивается осреднением по $N_b$ блокам и для данной частотной компоненты $f_k$ может быть представлена в виде $S_{f_k}={\widehat{Q}}_{f_k}^{}{\widehat{Q}}_{f_k}^{\ast }$, где ${\widehat{Q}}_{f_k}^{}=\sqrt{\Delta t/N_b}[q_k^{(1)}\mathrm{,...,}{q}_k^{(N_b)}]$ – матрица, состоящая из k-х столбцов матриц $Q^{(n)}$.

Поскольку матрица ${\tilde{S}}_{f_k}=S_{f_k}W$ является эрмитовой, для нее существует разложение по ортогональному базису ${\tilde{S}}_{f_k}=U_{f_k}{\Lambda }_{f_k}{U}_{f_k}^{\ast }$, которое можно выполнить стандартными алгоритмами, откуда следует

 ${\Psi }_{f_k}=W^{-1/2}{U}_{f_k}$,            (17)

что решает задачу SPOD-разложения. При этом найденные так собственные векторы ортогональны в смысле введенного скалярного произведения ${\Psi }_{f_k}^{\ast }W{\Psi }_{f_k}=I$ (матричная форма выражения (14)), а для взаимной спектральной плотности имеет место разложение

$S_{f_k}={\Psi }_{f_k}{\Lambda }_{f_k}{\Psi }_{f_k}^{\ast }$,            (18)

являющееся матричным аналогом выражения (15). В практическом смысле соотношение (18) означает, что интенсивность акустического излучения на заданной частоте $f_k$ (диагональные элементы матрицы $S_{f_k}$ ), полученная в дискретном наборе точек ${\theta }_i$, представляется в виде суперпозиции $N$ слагаемых

$S_{fk}={\displaystyle \sum _{j=1}^N{\lambda }_{fkj}}{\left|{\Psi }_{fkj}\right|}^2,$,     (19)

где $s_{f_k}$ – вектор, составленный из диагональных элементов матрицы $S_{f_k}$, ${\lambda }_{f_k}{}_j$ – найденные собственные числа – диагональные элементы матрицы ${\Lambda }_{f_k}$, ψ_fk j – найденные собственные векторы – столбцы матрицы ${\Psi }_{f_k}$. В соответствии с построенной процедурой старшие слагаемые в (19) определяют направленности и интенсивности базисных мультипольных источников. В качестве меры интенсивности каждого мультиполя можно принять величину в максимуме его излучения

w_fk j= λ_fk jmax _j [ψ^*_fk jψ_fk j].               (20)

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

Рассмотрим несколько модельных примеров с типичными параметрами лабораторных аэроакустических экспериментов. Пусть в начале координат задано три базисных диполя (т.е. степень мультиполя l = 1) с направленностями излучения, характеризующимися функциями $P_1^0(\cos \theta )=\cos \theta$, $P_1^1(\cos \theta )\cos \phi =\sin \theta \cos \phi$ и $P_1^1(\cos \theta )\sin \phi =\sin \theta \sin \phi$, с одинаковыми амплитудами, некоррелировано излучающих на частоте $f_0=1$ кГц (рис. 1). Некоррелированность в данном случае означает, что в каждой временной реализации конечной длительности начальные фазы каждого источника являются случайными величинами, равномерно распределенными на отрезке $[\mathrm{0,2}\pi )$, а математическое ожидание остальных случайных величин, представляющих интерес, вычисляется с помощью осреднения по реализациям случайного процесса. Рассмотрим точки наблюдения в дальнем поле, лежащие на дуге окружности радиуса R=1 м в плоскости, соответствующей азимутальному углу $\phi =0$ (рис. 1). Для указанных наблюдателей принимаемый сигнал в каждой временной реализации будет иметь вид

 $\begin{array}{l}{\widehat{p}}_{f_0}(R,\theta \mathrm{,0,}t)\approx \frac{e^{i2\pi f_0(t-R/c)}}{R}\times \\ \times (a_0{e}^{-i{\alpha }_0}\cos \theta +a_1{e}^{-i{\alpha }_1}\sin \theta ),\end{array}$  (21)

где $a_0$, $a_1$ – амплитуды источников, ${\alpha }_0$, ${\alpha }_1$ – фазы источников (различные в разных реализациях).

Рис. 1. Модельные источники из некоррелированных диполей (l =1) или квадруполей (l =2).

Аналогичный модельный источник может быть построен из квадруполей, излучение которых при $\phi =0$ будет описываться выражением

 $\begin{array}{c}{\widehat{p}}_{f_0}(R,\theta \mathrm{,0,}t)\approx \frac{e^{i2\pi f_0(t-R/c)}}{R}\left(a_0{e}^{-i{\alpha }_0}\frac{1}{2}(3{\cos }^2\theta -1)+\right.\\ \left.+a_1{e}^{-i{\alpha }_1}\frac{3}{2}\sin 2\theta +a_2{e}^{-i{\alpha }_2}3{\sin }^2\theta \right).\end{array}$ (22)

В модельных примерах с дипольными источниками рассмотрено три случая с различными соотношениями их амплитуд: $a_0$ = 2 и $a_1$ = 1, $a_0$ = 1 и $a_1$ = 2, $a_0$ = 1 и $a_1$ = 1. Отметим, что случай, когда амплитуды источников примерно равны, и суммарное их излучение практически изотропно, представляет наибольший интерес для анализа, поскольку в противном случае (явное доминирование определенного источника) характер мультиполя может быть определен непосредственно из формы диаграммы направленности излучения.

Сигналы задавались дискретно с шагом $\Delta t={10}^{-5}$ c, длина каждой реализации составляла $N_f=1024$, число реализаций $N_b=100$. Число точек наблюдения $N$ было выбрано равным от 10 до 30, что типично для аэроакустических экспериментов. Результаты работы SPOD-алгоритма приведены на рис. 2. На рис. 2а, 2б и 2в показаны результаты для различных амплитуд источников при $N=10$. Линиями представлены направленности исходно заданных источников и направленности суммарного шума, символами показаны SPOD-моды. Видно, что для всех случаев имеет место хорошее соответствие между SPOD-модами и исходными источниками как по амплитуде, так и по форме направленности. Число значимых SPOD-мод равно числу мультиполей данного порядка (остальные моды много меньше по амплитудам), старшая мода соответствует более интенсивному источнику. Увеличение числа наблюдателей приводит к улучшению соответствия (рис. 2г).

Рис. 2. Направленности излучения исходных дипольных источников (линии) и SPOD-мод (маркеры): 1 – осесимметричный диполь (n = 0); 2 – поперечный диполь (n = 1); 3 – суммарный шум. ○ – SPOD-мода j = 1; ◊ – SPOD-мода j = 2; □ – суммарная интенсивность SPOD-мод. (а) –  $a_0$= 2,  $a_1$= 1, N = 10; (б) – $a_0$ = 1,  $a_1$= 2, N = 10; (в) –  $a_0$= 1,  $a_1$= 1, N = 10; (г) – $a_0$ = 1,  $a_1$= 1, N = 20.

Для квадрупольных источников рассматривался случай, когда их суммарное излучение в среднем изотропно, что имеет место при $a_1=2a_2=a_0/\sqrt{3}$, a₀ = 2/3. Результаты SPOD-разложения представлены на рис. 3 для случаев N = 10 и N = 30. Как видно, при малом числе точек наблюдения амплитуды излучения каждой моды восстанавливаются достаточно точно, однако детали направленностей SPOD-мод в зонах, удаленных от их максимумов, могут отличаться от направленностей исходных источников. Это объясняется тем, что при увеличении степени мультиполей, а значит, числа лепестков в их направленностях, условие ортогональности в дискретной форме (14), записанное для фиксированного числа точек наблюдения N, будет все хуже аппроксимировать исходное условие (10). Это, в свою очередь, означает, что SPOD-моды, для которых, по построению, условие (14) выполняется строго, будут все более отличаться от направленностей исходных источников. Эту проблему решает уменьшение шага дискретизации по пространству (рис. 3б).

Рис. 3. Направленности излучения исходных квадрупольных источников (линии) и SPOD-мод (маркеры) для (а) – N = 10 и (б) – N = 30; 1 – осесимметричный квадруполь (n = 0); 2 – квадруполь (n = 1); 3 – квадруполь (n = 2); 4 – суммарный шум. ○ – SPOD-мода j = 1; ◊ – SPOD-мода j = 2; Δ – SPOD-мода j = 3; □ – суммарная интенсивность SPOD-мод.

2. Экспериментальная проверка метода

Для экспериментальной проверки разработанного алгоритма в заглушенной камере АК-2 ЦАГИ были проведены измерения шума обтекания цилиндра, при котором, как известно [1, 8–10], реализуются компактные источники дипольного типа, причем в шуме проявляются две дипольных компоненты, ортогональные оси цилиндра: так называемые диполи подъемной силы и силы сопротивления.

Стальной цилиндр диаметром d = 5 мм помещался в турбулентную струю, истекающую из круглого сопла диаметром D = 40 мм со скоростью V_j = 100 м/с. Цилиндр располагался в зоне турбулентного течения на расстоянии 250 мм вниз по потоку от среза сопла (рис. 4). Данная конфигурация соответствовала эксперименту, описанному в [33], что позволило провести сравнение результатов идентификации источников, полученных различными методами.

Рис. 4. Схема эксперимента. Стрелками показаны ориентации дипольных моментов шума обтекания цилиндра: диполя подъемной силы (L) и диполя силы сопротивления (D).

Шум измерялся с помощью 9 микрофонов Bruel&Kjaer (тип 4189), равномерно распределенных с шагом Δθ = 15° по дуге окружности радиуса R = 1.08 м с центром в точке пересечения цилиндра с осью струи, как показано на рис. 4. Измерения проводились как для свободной струи, так и для струи с цилиндром. Цилиндр ориентировался так, чтобы его ось была ортогональна плоскости расположения микрофонов (рис. 4). В таком случае дипольные моменты, связанные с пульсациями подъемной силы и силы сопротивления (показаны стрелками на рис. 4), лежат в плоскости микрофонов, и на массиве микрофонов реализуется звуковое поле, близкое к полю, описываемому выражением (21).

При наличии цилиндра в общем шуме такой системы доминирует именно шум цилиндра, как показано на рис. 5, поскольку шум струи носит квадрупольный характер, и в низкоскоростных течениях его интенсивность мала по сравнению с интенсивностью дипольных источников [1]. На спектрах на рис. 5 число Струхаля рассчитано по диаметру цилиндра и скорости струи $\text{St}=fd/V_j$. Максимум шума цилиндра соответствует числу Струхаля около St = 0.14, что соответствует результатам измерений работы [33].

Рис. 5. Измеренные спектры шума: (а) – θ=30°; (б) – θ=90°; (в) – θ=150°. 1 – цилиндр в струе; 2 – струя без цилиндра.

На рис. 6 приведены спектры собственных чисел SPOD-разложения звукового поля при наличии цилиндра в струе (моды упорядочены по убыванию энергии, число мод равно числу микрофонов). Как видно, практически во всем частотном диапазоне доминируют первые две SPOD-моды, что косвенно указывает на дипольный характер шума. Интересно отметить, что при St ≈ 0.06 интенсивности двух мод сравниваются, и спектр первой моды практически непрерывно переходит в спектр второй моды и наоборот. Данный эффект также указывает на то, что SPOD-моды, определенные по предложенному в настоящей работе алгоритму, соответствуют физическим объектам, связанным с определенными механизмами шумообразования.

Рис. 6. Спектры собственных чисел SPOD-разложения звукового поля при наличии цилиндра в струе.

Связь SPOD-мод с физическими мультиполями подтверждается и структурой их направленности, изображенной на рис. 7 для различных полос частот. Показаны только две доминирующие моды. Видно, что одна из них имеет однолепестковую структуру, другая – двухлепестковую (причем соседние лепестки всегда находятся в противофазе), что соответствует различным компонентам дипольного излучения, доминирование которого ожидается в шуме обтекания цилиндра. В низкочастотной области St < 0.06 доминирует излучение продольного диполя, связанного с пульсациями силы сопротивления (рис. 7а), а на частотах выше St = 0.06, в том числе, в области спектрального пика, определяющим становится поперечный диполь, соответствующий пульсациям подъемной силы (рис. 7б-7г). Отклонения форм направленностей от “идеальных” диполей могут быть связаны с эффектами конвекции и рефракции, “деформирующими” направленности исходных базисных источников, что может приводить к некоторому сдвигу положений минимумов направленностей в силу условия взаимной ортогональности SPOD-мод (рис. 7а), а также с малостью амплитуды одного из источников, находящейся на границе динамического диапазона, в результате чего характеристики его направленности восстанавливаются с заметной погрешностью (рис. 7в, 7г). При этом отметим, что конвективное усиление суммарного шума вниз по потоку, хорошо заметное в общем шуме (рис. 7а), в целом достаточно корректно переносится алгоритмом SPOD на соответствующие дипольные компоненты, по крайней мере, в зоне максимумов направленности.

Рис. 7. Направленности излучения в полосах частот для цилиндра в струе: (а) – St = 0.017; (б) – 0.08; (в) – 0.14; (г) – 0.2. Сплошная линия – измерения; ○ – SPOD-мода j = 1; ◊ – SPOD-мода j = 2; □ – суммарная интенсивность SPOD-мод.

Количественной мерой излучения каждой компоненты, как отмечалось выше, можно считать ее интенсивность в максимуме излучения (для SPOD-мод это определение соответствует формуле (20)). Аналогичная величина была выбрана в работе [33], что позволяет провести прямое сравнение выделения спектральных характеристик интенсивности дипольных источников, полученных в [33] методом азимутальной декомпозиции, и методом SPOD-разложения, предложенным в настоящей работе. Такое сравнение приведено на рис. 8, и оно демонстрирует совпадение результатов оценки интенсивностей дипольных компонент, полученных различными методами. При этом метод на основе SPOD-разложения использует одномерный массив микрофонов, расположенный в одной азимутальной плоскости, что в некоторых случаях делает его проще реализуемым по сравнению c методом азимутальной декомпозиции или методом плоского бимформинга [33].

Рис. 8. Спектральные плотности интенсивности излучения отдельных SPOD-мод (1-9), вычисленные по формуле (20), и спектральные плотности излучения поперечного (10) и продольного (11) диполей, определенные с помощью метода азимутальной декомпозиции [33] (данные приведены к расстоянию 1 м от источника).

На рис. 9 приведены результаты SPOD-разложения для звукового поля свободной струи. Отметим, что для анализа шума струи используется, как это обычно принято, число Струхаля ${\text{St}}_j=fD/V_j$, рассчитанное по скорости струи и диаметру сопла. Данная струя является низкоскоростной, и максимум спектра ее излучения приходится на частотный диапазон 300–600 Гц, что соответствует St_j = 0.1…0.2. В области максимального излучения, как видно из спектров собственных чисел SPOD-разложения, показанных на рис. 9а, в звуковое поле вносят вклад три моды, что указывает на квадрупольный характер излучения [11]. Отметим, что предложенный в настоящей работе вариант SPOD-метода предполагает, что источник звука является компактным и расположен в центре дуги, на которой находятся микрофоны. При этом, как известно [3], характерное положение источников шума в турбулентных струях зависит от частоты их излучения так, что более высокочастотные источники располагаются ближе к соплу. Поэтому для использованной в эксперименте дуги микрофонов, центрированной на точке, удаленной от среза сопла на расстояние $6.25D$, пространственные формы отдельных SPOD-мод будут близки к физическим мультиполям для той частоты, которой соответствует расположение источников в области центра микрофонного массива. Для указанного расстояния от среза сопла характерные частоты излучения соответствуют как раз максимуму шума струи, т.е. частотному диапазону St_j = 0.1…0.2 [3, 24]. При этом, поскольку источники шума струи даже для узкого диапазона частот не являются компактными в продольном направлении [24], соответствие направленностей SPOD-мод направленностям точечных квадруполей будет не вполне точным. Это видно из рис. 9б, на котором приведены направленности излучения доминирующих SPOD-мод для St_j = 0.1. Видно, что качественно формы SPOD-мод соответствуют квадруполям различного азимутального порядка, формы которых могут быть с высокой точностью получены с помощью метода азимутальной декомпозиции [11, 27], с некоторым отличием в расположении минимумов, связанным с конвективно-рефракционным усилением излучения вниз по потоку [20] (обсуждавшимся выше применительно к шуму цилиндра), а также в величине амплитуд тех лепестков, которые находятся на границе динамического диапазона метода (~10 дБ ниже доминирующей моды), что отчасти связано с некомпактностью источников. Тем не менее, общая структура SPOD-мод позволяет вполне определенно отнести их к соответствующим квадрупольным источникам, т.е. адекватно оценить доминирующий механизм генерации шума исследуемым течением.

Рис. 9. (а) – Cпектры собственных чисел SPOD-разложения звукового поля свободной струи; (б) – направленности излучения доминирующих SPOD-мод для St_j = 0.1; обозначения кривых как на рис. 3.

Заключение

В данной работе предложено использовать метод собственного ортогонального разложения в частотной области (SPOD) для идентификации мультипольной структуры аэроакустических источников по измерениям в дальнем звуковом поле. Преимуществом SPOD-метода является то, что он может быть реализован в рамках линейного микрофонного массива, располагаемого в одной азимутальной плоскости, в то время как метод азимутальной декомпозиции, с помощью которого такая задача может быть решена наиболее точно, требует использования азимутальных круговых микрофонных массивов, заметающих существенную долю телесного угла, в который излучается шум. Показано, что при условии компактности источника определенного (заранее неизвестного) мультипольного порядка, и расположении его в центре кругового микрофонного массива, заметающего сектор, близкий к 180°, можно построить такую процедуру SPOD-разложения за счет выбора скалярного произведения, в которой доминирующие SPOD-моды соответствуют физическим мультипольным компонентам излучения.

Метод протестирован на примерах с модельными точечными мультиполями, а также проверен в эксперименте применительно к анализу шума обтекания цилиндра и шума струи. Эффекты конвекции и рефракции снижают точность восстановления направленностей отдельных мультиполей, однако при малых числах Маха потока корректность определения интенсивностей доминирующих мультипольных компонент сохраняется. При наличии теоретических моделей влияния конвекции и рефракции на диаграмму направленности излучения базисных мультиполей скалярное произведение в методе SPOD можно скорректировать, повысив тем самым точность восстановления мультипольной структуры.

Отметим в заключение отличие разработанного метода на базе SPOD от метода бимформинга, использованного, например, в работах [31–33] для идентификации дипольных источников шума. Оба подхода связаны с анализом матрицы взаимных спектров, измеряемой некоторым массивом микрофонов. В методах типа бимформинга такая матрица моделируется с помощью исходно выбранных моделей источников заданной мультипольности. В таком подходе корректность результата зависит от адекватности набора источников, заложенных в алгоритм. Преимуществом его является возможность определения заранее неизвестного местоположения источников шума, недостатком – необходимость априори определять тип источников. В предложенном методе типа SPOD матрица взаимных спектров раскладывается по эмпирическому ортогональному базису, что всегда оказывается возможным, и при определенных условиях компоненты этого базиса могут соответствовать физически осмысленным мультиполям. Преимуществом данного метода является возможность определения порядка мультипольности источника в результате измерений, недостатком – необходимость априорного знания положения источника. Ограничением SPOD-метода является также требование компактности источников шума: чем хуже выполняется условие компактности, тем сложнее интерпретируются результаты.

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

Экспериментальные данные, использованные в работе, получены в УНУ “Заглушенная камера с потоком АК-2” ФАУ “ЦАГИ”, модернизированной при поддержке Министерства науки и высшего образования Российской Федерации.

Работа выполнена при финансовой поддержке РНФ (грант 19-71-10064).