Применение акустической аналогии для расчета звука в дальнем поле

Приведем вывод волнового уравнения Фокс Вильямса–Хоукингса, следуя работе [1], для общего случая, когда поверхность, внутри которой распределены источники (поверхность Кирхгоффа), предполагается подвижной. Пусть f ( X,t ) = 0 представляет собой уравнение подвижной поверхности, охватывающей регион, в котором происходит генерация звука (рис.1), x ^r – пространственная координата, t – время.

Рассмотрим задачу нахождения звукового поля в открытой области внешней к поверхности f (x, t) = 0 . Примем, что f > 0 во внешней области и f < 0 внутри поверхности, тогда функция Хэвисайда H( f ) = 1, если f > 0 и H( f) = 0, если f < 0 . Кроме того, предположим, что f (X, t) определена так, что Vf = n - внешняя нормаль к поверхности f (X, t)= 0. Это предположение не является ограничением общности, так как всегда можно переопределить f как f ^Vf. Следует отметить, что поверхность f (X, t) = 0 не предполагается непроницаемой. Пусть, кроме того, поверхность может двигаться со скоростью с компонентами vi так, что af(.r, t)+v af(x, t) = 0

a t i a x i         '

В этом случае

a H(f) _a h a f = _ а н a f a t "af a t ” ^v a f a x,

а н

d f

= ^_ v ^ ' f )

d x-

n = V f

f >  0

. f <  0

f(x,t) = 0

( 1

( c 2 d t²

d 2

V ² 1 p ' = 0

Рис.1. Подвижная деформируемая поверхность f ( x , t ) = 0

Продолжим все поля газодинамических параметров внутри поверхности f ( x, t ) = 0 величинами, относящимися к невозмущенному потоку. Это продолжение задачи на все неограниченное пространство необходимо, чтобы при решении волнового уравнения можно было использовать функцию Грина волнового уравнения для неограниченного пространства, которая имеет простой вид. Таким образом, при прохождении через поверхность f ( x, t ) = 0 все газодинамические параметры будут терпеть разрыв. Предположим, кроме того, что параметры потока не имеют никаких других разрывов, кроме разрывов при переходе через поверхность f ( x, t ) = 0.

Уравнение непрерывности и сохранения импульса во внешней области f > 0 при отсутствии массовых сил имеют вид

Sp+ d(pu1) = 0 , dt

d(PUi) + d(pUUj + P8ij -Tij 1 = 0 dtd

где	p - плотность,		p	– давление,	ui   – компоненты   скорости,
L , =И	<5 u d u , —L + —- (a x j a x i 7	2 . s uk --8 ij------ 3    d x k _	-	тензор вязких	напряжений, p - коэффициент

динамической вязкости, 8 j - дельта Кронекера.

Умножив уравнения (3) и (4) на функцию Хэвисайда и используя выражение (2), получим

((p - p0 )H(f)) +    (pUH(f)) = p0vi jH + p(ui - vi) IH ,(5)

dt                     dxi                      dxid

I" ⁽ p u i H ^{( f} )) ^{+     ((} p uU j ⁺ P ⁸ ij ^{- T} ij ) H ^{( f} )) = ⁽ p ui ⁽ u j ^- v j ) ⁺ p ⁸ ij ^{- T} ij ) | H , (6)

dt               dxid где p0 - средняя плотность.

Уравнение Фокс Вильямса–Хоукингса может быть получено из (5) и (6) следующим образом. Сначала возьмем производные по времени от обеих частей уравнения (5), затем возьмем дивергенцию от обеих частей уравнения (6) и потом вычтем последнее уравнение из первого. Наконец, вычтем V2(c2(p-p0)H(f)) из обеих сторон полученного уравнения и перегруппируем полученное выражение так, чтобы волновой оператор, действующий на c2 (p - p0 )H(f), стоял в левой части. Таким образом, получим уравнение, которое называют уравнением Фокс Вильямса– Хоукингса для подвижной проницаемой поверхности:

' 1 d ²

c 2 a/² V ^c 0 ^d t

^

V ² ( c 2 ( p-p 0 ) H ( f )) =

——(TijH(f ))+ — ((p0Vn +p(Un - Vn))8(f ))- dxidxj              d t                                    (7)

^- d x - (( ( P ^- P 0 ) n i ^-T ij n j ⁺p ^ui ⁽ u n ^{- V}n ^{)) 8 ( f} ) )

Здесь введены следующие обозначения: n i = f /d x i - компоненты внешней единичной нормали к f = 0 u_n = u i n i - локальная скорость среды в направлении нормали, v n = v i n i - локальная нормальная скорость поверхности f = 0, T ij = p u i u j + ( ( p - p 0 ) - c 2 ( p - p 0 ) ) s ij - t ij - тензор Лайтхилла, c ₀ - скорость звука.

Из уравнения (7) видно, что правая часть зависит от неизвестных скоростей и давлений в среде, поэтому, строго говоря, не может рассматриваться как внешний заданный источник и должна быть определена из общих уравнений аэродинамики. Лайтхиллом [2] предложено рассматривать правую часть уравнения (7) как заданное эквивалентное нестационарному потоку распределение акустических источников, которые излучают в идеальную неподвижную среду. Это допущение получило название акустической аналогии Лайтхилла. Его применение возможно, когда несущественно взаимодействие звуковых волн с потоком.

Будем считать, что расположение поверхности Кирхгоффа выбрано так, что с ее наружной части распространение возмущений происходит линейно, т.е. (p — p0) = c2(p —Po), кроме того, и сама поверхность находится в линейной зоне, поэтому вкладом тензора вязких напряжений и тензора Лайтхилла можно пренебречь.

Ограничимся теперь рассмотрением неподвижной проницаемой поверхности f = 0 . С учетом этих предположений уравнение (7) примет вид

f ,

V С 02 d t ²

х

V ² p’H ( f) = -|(( p U n ) S ( f ) ) d t

d

^— —^{(( (} p ^— p 0 ) nt + P uiu n ) 3 ⁽ f ) ) ,     (8)

a x i

где через p' = p — p 0 обозначено звуковое давление.

Решение уравнения (8) может быть найдено аналитически с использованием функции Грина для свободного пространства [3].

Рассмотрим волновое уравнение следующего общего вида

f 7° V c 0 a t²

V ² ] p ' = Q ( y, т ).

Решение уравнения (9) может быть получено в виде

TO p ‘ (x, t) = J J Q(y, t) G(x, 11 y, T)d3 ydт ,

T —TO где G(x, t1 y, t) - функция Грина, которая удовлетворяет волновому уравнению для импульсного точечного источника, расположенного в точке x = y и излучающего импульс в момент времени t = т , т.е.

1 а ²

( c 0 ² a t ²

^{— V} 2

G = 5 ( x — y ) 5 ( t — т ).

Кроме того, на функцию G накладывается следующее условие: G = 0 при t < т .

Для свободного пространства существует единственная

удовлетворяющая уравнению (11) и условию (12). Она имеет вид

G ( x , 1 1 y , т ) =

4л| x — y|

f

^{8 t}

V

— т —

^c 0 7

(12) функция,

Используя (13), можно получить следующий вид решения уравнения (8):

4np' (x, t) = J f=0 -

P^ dS + ±J L. dS +J r J      c 0 f=0 - r J      f=0-

IdS

,

7 x

где введены обозначения U _n

p u i

--- n ,

V ^p 0 7

L r =^{( (} p ^— p 0 ) n i ⁺ P uiu n ) r i , r i

rr

x — y

M "

компоненты единичного вектора в направлении от точки на поверхности Кирхгоффа до точки наблюдения, r = Jr — y - расстояние между точкой наблюдения и точкой на поверхности Кирхгоффа. Точка над величинами в (14) означает производную по времени, а величины в квадратных скобках берутся в момент времени т = t — r]cq .

Для нахождения интеграла (14) используются методы интегрирования. Рассмотрим интеграл следующего общего вида:

4np' (x, t) =

I f=0 l

Q ( У , т )

r

dS

,

численного

где Q(y, т ) - заданная функция положения источника и времени.

Интегрирование в формуле (15) происходит по поверхности f = 0, а интегранд рассматривается во время излучения звука, т.е. в во время т = t — r/c q .

Наиболее простым и распространенным методом численного определения интеграла (15) является следующий. Поверхность f = 0 разбивается на набор плоских элементов, для каждого из которых интегранд вычисляется в центральной точке элемента, в соответствующее этой точке запаздывающее время и умножается на площадь элемента [4]. После этого производится суммирование по всем элементам. Таким образом, исходный интеграл аппроксимируется выражением

4пp' (x, t) «f i=1

Q У, , t — V

где N – общее количество элементов, y ^r _i – центральная точка i -й панели.

Тестовые задачи

В качестве первой тестовой задачи для верификации правильности алгоритма была использована задача об излучении сферических волн точечным источником.

Точечный источник

Точка наблюдения

100 м

R=1 м

Рис.2. Задача о точечном источнике

Аналитическое решение для этой задачи в сферических координатах имеет вид p ' =

^^^^^^^B

p 0 ю A

r

^^^^^^^B

t

u =

^^^^^^^^

ю A . r

--Sin ю — c 0 r V V c 0

В качестве поверхности Кирхгоффа была взята сфера радиусом 1 метр (рис.2), на которой в соответствии с формулами (17) задавались пульсации давления плотности и скорости. Расчеты были выполнены для двух вариантов разбиения сферы на элементы. Для первого варианта расчетная сетка составляла 600 (6x10x10) элементов, для второго – 5400 (6x30x30). Точка наблюдения находилась на расстоянии 100 метров от центра сферы. Значения констант, при которых проводились расчеты, были взяты следующие: A =10 Па, ω=2000π рад/с, p0 =101325 Па, ρ0 =1,22 кг/м3, c0=340 м/с. Полученные сигналы в сравнении с аналитическим решением представлены на рис.3. Из рисунка видно, что при уменьшении размеров элементов решение приближается к аналитическому и для второго варианта практически с ним совпадает. Отношение длины волны излучаемого звука к линейному размеру элемента для второго варианта составляет примерно 6,5. Таким образом, можно сделать оценку для выбора густоты сетки на поверхности Кирхгоффа: для правильного описания дальнего звукового поля необходимо около семи элементов на длину волны излучения.

Рис.3. Пульсации давления на расстоянии 100 метров от источника

В качестве второй тестовой задачи для применения алгоритма (16) численного интегрирования выражения (14) был выбран случай излучения звука при обтекании цилиндра потоком вязкого газа. Условия проведения эксперимента и спектр экспериментально измеренного уровня звукового давления были взяты из статьи [5]. Для нахождения нестационарных полей давления, скорости и плотности вокруг цилиндра был использован газодинамический пакет Fluent версии 6.1. Вид расчетной области показан на рис.4. Количество узлов расчетной сетки составило 303000.

При проведении газодинамического расчета задавались следующие граничные условия: на входе (рис.4) – скорость набегающего потока 15 м/с, температура 290 К; на выходе и верхней и нижней границах – открытая граница с атмосферным давлением (101325 Па) и температурой 290 К. Диаметр цилиндра – 10 мм, длина цилиндра в расчете взята равной трем диаметрам (реальная длина в эксперименте 50 см.). На боковых поверхностях ставились условия периодичности. Число Рейнольдса для заданного режима течения составляет порядка 10 ⁴ .

При расчете использовался неявный метод контрольных объемов со вторым порядком точности по пространству и по времени и LES подход для моделирования турбулентности. Газодинамический расчет проводился в два этапа, сначала решение выводилось на периодический режим. Для примера на рис.4 приведено полученное в расчете мгновенное поле скорости в среднем по длине цилиндра сечении. Затем производилась запись в файлы нестационарных полей давления, плотности и скорости на четырех выбранных поверхностях с шагом по времени 0,0001 с в течение 1000 шагов. Выбранные поверхности Кирхгоффа представляют собой цилиндрические поверхности с радиусом, равным 0,5; 1, 1,5 и 2 диаметра цилиндра, и осью, совпадающей с осью цилиндра (см. рис.4). Будем называть их соответственно 0,5D, 1D, 1,5D и 2D.

Рис.4. Вид расчетной области

Рис.5. Звуковое давления в точке наблюдения

Нестационарные поля на поверхностях 0,5D, 1D, 1,5D и 2D последовательно использовались для численного вычисления интеграла (14) в точке, где производился замер уровня звука в эксперименте. Эта точка находится на расстоянии 1 м над цилиндром, в средней точке по его длине, в плоскости, проходящей через ось цилиндра. Сигнал в этой точке, полученный для различных поверхностей Кирхгоффа, представлен на рис.5.

Из рисунка видно, что чем больше радиус, на котором располагается поверхность Кирхгоффа, тем больше гармоник присутствует в сигнале. Это связано с тем, что при увеличении радиуса все больше источников охватывается поверхностью Кирхгоффа. Начинают учитываться не только дипольные источники, связанные с силовым взаимодействием цилиндра и потока (как в случае совпадения поверхности Кирхгоффа с поверхностью цилиндра 0,5 D ), но также и квадрупольные источники, связанные с турбулентностью в охваченной поверхностью Кирхгоффа области течения.

Для полученных сигналов было выполнено преобразование Фурье и найден уровень звукового давления в децибелах. Результат представлен на рис.6. Из рисунка видно, что результаты расчета качественно правильно воспроизводят основную частоту срыва вихрей. Экспериментальное выражение для основной частоты срыва вихрей имеет вид [6]

и f = ShD, (18) где Sh * 0.2 - число Струхаля для цилиндра, и - скорость потока, D - диаметр цилиндра. Используя (18), получим значение основной частоты f ~ 300 Гц, что совпадает с основной частотой в эксперименте и хорошо согласуется с полученной численно (f ~ 310 Гц).

Частота, Гц

Рис.6. Уровень звукового давления в точке наблюдения

Расхождение по величине уровня звукового давления на основной частоте можно объяснить тем, что в расчетах, с целью экономии вычислительных ресурсов длина цилиндра была взята значительно меньше, чем в реальном эксперименте, следовательно, суммарная мощность источников получилась меньше. Для выяснения этого вопроса воспользуемся теорией Филлипса для оценки интенсивности звука при обтекании цилиндра [3]. Согласно ей для “коротких” цилиндров, т.е. для которых можно считать, что по всей длине цилиндра срыв вихрей происходит синфазно, имеет место следующее выражение для интенсивности звука в точке, лежащей в плоскости, проходящей через ось цилиндра перпендикулярно течению газа на расстоянии r от оси цилиндра,

I = к ² Sh ² l ² p 0 u ⁶

1       32c03r2     , где l - длина цилиндра, к - константа, принимающая значения от 0,5 до 1. Принимая Sh =0,2, p0=1,22 кг/м3, u = 15 м/с, c0 =340 м/с, r=1 м для значений к =0,5 и к = 1, получим следующие оценки интенсивности звука для цилиндра, использовавшегося в расчете (l =0,03м): I1к=0,5=0,99848е-7 Вт/м2, I1к=1=0,39939е-6 Вт/м2. Для сравнения с расчетом пересчитаем интенсивности в уровни звукового давления. Согласно формуле L = 10log(I/10), где I0=1е-12 Вт/м2 - предел слышимости, получаем уровни: LK=0,5=5O дБ, LK=0,5=56 дБ, то есть примерно такие же, как получаются в расчете Lрасч =54 дБ.

Сравним теперь, как согласуется теория Филлипса с экспериментальными данными. Для “длинных” цилиндров, т.е. для которых длина цилиндра велика по сравнению с расстоянием, на котором срывающиеся с цилиндра вихри коррелированы,

I = 17 D к ² Sh ² l p ₀ u ⁶

.

• ²         32 c ₀ ³ r ²

Вычислим интенсивность звука для цилиндра, использовавшегося в эксперименте (l =0,5м). Для к =0,5 и к = 1 получим 12к=0,5 =0,943e-5 Вт/м2, I2х 1 0,377е4 Вт/м2, что соответствует уровням звука Lк=0,5 =69,75 дБ и Lк=1=75,77 дБ. Экспериментальный уровень на основной частоте срыва вихрей составляет Lэксп =68 дБ (см. рис.6), т.е. теория несколько завышает уровень звукового давления, но все равно достаточно хорошо согласуется с экспериментом.

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

Рис.7. Направленность на частоте 310 Гц

Сравнение результатов, полученных для разных поверхностей Кирхгоффа, показывает, что для всех поверхностей уровень звукового давления на основной частоте срыва вихрей получается примерно одинаковый. Различия проявляются в предсказании гармоник основной частоты. Наиболее близки к экспериментальным данным результаты, полученные с использованием поверхности 1 D . Для этой поверхности была проведена серия расчетов для нахождения направленности излучения. Точки наблюдения располагались на радиусе 1 м через 10 градусов в плоскости, перпендикулярной оси цилиндра, проходящей через его середину. Зависимость уровня звукового давления от угла для основной частоты срыва вихрей (310 Гц) представлена на рис.7.

Как видно из графика, направленность имеет форму восьмерки, что характерно для дипольных источников звука.