1. На главную
  2. Электронные версии
  3. Акустический журнал
  4. T. 69, № 6, 2023

Акустический журнал, 2023, T. 69, № 6, стр. 713-721

Метод конечно-элементного моделирования гидродинамического шума, возникающего при обтекании упругих тел

А. С. Суворов a, Е. М. Соков a, А. Л. Вировлянский a, В. О. Еремеев a, *, Н. В. Балакирева a

a Институт прикладной физики РАН
603950 Нижний Новгород, ул. Ульянова 46, Россия

* E-mail: eremeev.vladimir.o@ipfran.ru

Поступила в редакцию 17.04.2023
После доработки 16.08.2023
Принята к публикации 19.09.2023

Полный текст (PDF)

Аннотация

Представлен конечно-элементный метод расчета гидродинамического шума, возбуждаемого турбулентными пульсациями жидкости в присутствии упругого тела. Традиционный подход к решению этой задачи на основе прямого решения уравнения Лайтхилла требует большого объема вычислений. Показано, что ситуация существенно упрощается при расчете компонент шума на относительно низких частотах, которые отвечают длинам волн, превышающим размеры турбулентной области. В этом случае шумовое поле удается выразить через давление турбулентных пульсаций на поверхности упругого тела, найденное в приближении несжимаемой жидкости. Статья подготовлена по материалам доклада, представленного на IX Российской конференции “Вычислительный эксперимент в аэроакустике и аэродинамике”, г. Светлогорск, 26 сентября–1 октября 2022 г.

Ключевые слова: гидродинамический шум, конечно-элементное моделирование, излучение звука

ВВЕДЕНИЕ

В данной работе рассматривается одна из классических задач аэро- и гидроакустики: расчет звукового поля, возбуждаемого турбулентными пульсациями скорости с малыми числами Маха в присутствии упругого тела. Данный класс задач успешно решается методом физического эксперимента, однако он отличается дороговизной и требует много времени и подготовки к проведению эксперимента. Некоторые сложности, связанные с проведением физического эксперимента, изложены в [1–3]. В связи с этим большую популярность набирают численные эксперименты, что позволяет в настоящее время решать задачи, связанные с анализом генерации звука турбулентностью, с помощью так называемого гибридного подхода [4, 5]. Он базируется на предположении о том, что возбуждаемое звуковое поле столь слабо, что его обратным влиянием на турбулентность можно пренебречь. Задача решается в два этапа.

На первом этапе турбулентные пульсации рассчитываются путем численного решения уравнения Навье–Стокса. В зависимости от поставленной задачи и располагаемых вычисленных ресурсов используют либо метод моделирования крупных вихрей (Large Eddy Simulation, LES) [6–8], позволяющий описать вихревые структуры с минимальным размером, связанным с размером ячейки сеточной модели, либо RANS, описывающий только крупные энергонесущие вихревые структуры.

На втором этапе результаты этих вычислений применяются для оценки генерируемых звуковых полей. Это делается с помощью акустической аналогии Лайтхилла [9], подхода Фокса Вильямса и Хоукингса [10], соотношений, описывающих генерацию звука вихрями (vortex sound) [11] и многочисленных методов, основанных на использовании гидродинамической теории возмущений [3, 12]. Пример применения одного из гибридных подходов приведен в [13].

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

Учет деформации упругого тела требует применения другого подхода. В работах [15, 16] описано решение данной задачи с использованием метода конечных элементов. На этапе гидродинамического расчета определяются источники гидродинамического шума в виде тензора Лайтхилла, при этом поверхность обтекаемого тела считается абсолютно твердой. При расчете генерируемого звукового поля, на основе определенных источников, предполагается, что давление в жидкости описывается уравнением Лайтхилла, а смещения частиц среды внутри упругого тела задаются обычными уравнениями линейной теории упругости [17]. Совместное решение этих уравнений выполняется методом конечных элементов [15, 16]. Преимуществом данного метода, кроме учета вклада упругости конструкции в акустическое поле, является возможность подавления влияния граничных условий гидродинамической части задачи за счет наложения плавного пространственного фильтра. Кроме того, в классе задач, где относительно локальная неоднородность обтекаемого тела является интенсивным источником шумоизлучения, применяют подход с искусственной турбулизацией потока на входной границе расчетной области. В таких задачах без применения пространственной фильтрации искусственных турбулентных источников данные источники приводят к значительному искажению определяемого акустического поля. Однако практическая реализация данного подхода требует большого объема вычислений, так как расчет акустического поля базируется на использовании массива данных, включающих значения поля скорости во всех узлах гидродинамической модели среды для всех временных отсчетов.

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

ПОСТАНОВКА ЗАДАЧИ

Рассматривается ситуация, когда упругое тело – объем Ve, заполненный однородной упругой средой, – находится внутри турбулентной жидкости. Жидкость занимает все бесконечное пространство V вне области Ve. Турбулентные пульсации скорости жидкости заданы полем v(x, t), где x – радиус-вектор точки среды, t – время. Декартовы координаты векторов v и x обозначим соответственно vi и xi, i = 1, 2, 3. Компоненты вектора деформации в области Ve обозначим ui(x,t). Границу раздела жидкой и упругой сред обозначим символом S.

Предполагается, что турбулентные пульсации сосредоточены в некоторой окрестности упругого тела, характерный размер которой обозначим D. Нашей целью является расчет возбуждаемого турбулентностью гидродинамического шума на дистанциях $r \gg D$ от центра турбулентной области. В соответствии со сказанным во Введении, решение этой задачи для турбулентности с малыми числами Маха выполняется в два этапа.

1. Гидродинамический этап. Турбулентное движение жидкости численно рассчитывается в предположениях о том, что жидкость несжимаема, а поверхность упругого тела S является твердой. Результатом расчета являются функции vi(x, t) и pf(x, t), представляющие компоненты скорости и давления в турбулентной жидкости. Далее при выводе аналитических соотношений мы учитываем то обстоятельство, что нормальная производная давления pf на твердой поверхности S пропорциональна Re–1, где Re – число Рейнольдса [15, 18]. Ограничиваясь рассмотрением ситуаций, когда числа Рейнольдса достаточно велики, мы, следуя [18], приближенно считаем нормальную производную давления на поверхности S равной нулю.

2. Акустический этап. Расчет поля звукового давления p(x, t) путем решения уравнения Лайтхилла

(1)
$\frac{1}{{{{c}^{2}}}}\frac{{{{\partial }^{2}}p}}{{\partial {{t}^{2}}}} - \Delta p = \frac{{{{\partial }^{2}}{{T}_{{ij}}}}}{{\partial {{x}_{i}}\partial {{x}_{j}}}},$
где Tij(x, t) – тензор Лайтхилла, который в случае малых чисел Маха выражается приближенным соотношением [1, 2]
(2)
${{T}_{{ij}}} = {{\rho }_{0}}{{v}_{i}}{{v}_{j}},$
где ρ0 – плотность невозмущенной жидкости. Здесь и далее по повторяющимся индексам выполняется суммирование. Оставаясь в рамках традиционного подхода [4], правую часть уравнения Лайтхилла (1) рассматриваем как заданный источник звукового поля. Она находится подстановкой в (2) поля скорости турбулентных пульсаций, найденного на этапе 1.

Уравнение (1) для давления p(x, t) следует решать совместно с линейным уравнением для компонент вектора деформации (смещения) частиц упругого тела ui(x, t), i = 1, 2, 3 [17]. Связь между p и ui  устанавливается приведенными ниже граничными условиями на поверхности S.

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

$\tilde {p}\left( {{\mathbf{x}},\omega } \right) = \int\limits_{ - \infty }^\infty {p\left( {x,t} \right)} {\kern 1pt} {{e}^{{i\omega t}}}dt.$

Линейное волновое уравнение (1) для давления p переходит в уравнение Гельмгольца

(3)
${{\Delta }}\tilde {p}\left( {\mathbf{x}} \right) + {{k}^{2}}\tilde {p}\left( {\mathbf{x}} \right) = - \frac{{{{\partial }^{2}}{{{\tilde {T}}}_{{ij}}}\left( {\mathbf{x}} \right)}}{{\partial {{x}_{i}}\partial {{x}_{j}}}},$
где k = ω/c, c – скорость звука в жидкости. Здесь и далее аргумент ω для краткости опускаем. Внутри упругого тела уравнения для компонент вектора смещения [17] имеют вид
(4)
$ - {{{{\omega }}}^{2}}{{{{\rho }}}_{e}}{{\tilde {u}}_{i}}\left( {\mathbf{x}} \right) = \frac{{\partial {{{\tilde {\sigma }}}_{{ij}}}\left( {\mathbf{x}} \right)}}{{\partial {{x}_{j}}}},$
где ρe – плотность упругой среды, ${{{{\tilde {\sigma }}}}_{{ij}}}$ – компоненты тензора напряжений. В самом общем виде компоненты тензора напряжений в линейном приближении выражаются через компоненты вектора деформации соотношением [17]
(5)
${{{{\tilde {\sigma }}}}_{{ij}}} = {{A}_{{iji{\kern 1pt} 'j{\kern 1pt} '}}}\frac{{\partial {{{\tilde {u}}}_{{i'}}}}}{{\partial {{x}_{{j{\kern 1pt} '}}}}},$
где ${{A}_{{iji{\kern 1pt} 'j{\kern 1pt} '}}}$ – постоянные коэффициенты.

Уравнения (3) и (4) следует решать совместно с граничными условиями, заданными на поверхности S. Эти условия выражают непрерывность нормальной компоненты скорости

(6)
${{\left[ {\frac{{\partial{ \tilde {p}}\left( {\mathbf{x}} \right)}}{{\partial {{x}_{i}}}}{{n}_{i}}\left( {\mathbf{x}} \right) - {{{{\omega }}}^{2}}{{\rho }}{{{\tilde {u}}}_{i}}\left( {\mathbf{x}} \right){{n}_{i}}\left( {\mathbf{x}} \right)} \right]}_{{{\mathbf{x}} \in S}}} = 0$
и нормальной компоненты тензора напряжений

(7)
${{\left[ {{{{{{\tilde {\sigma }}}}}_{{ij}}}\left( {\mathbf{x}} \right){{n}_{j}}\left( {\mathbf{x}} \right) + \tilde {p}\left( {\mathbf{x}} \right){{n}_{i}}\left( {\mathbf{x}} \right)} \right]}_{{{\mathbf{x}} \in S}}} = 0.$

Здесь ni(x), i = 1, 2, 3 – компоненты единичного вектора n(x), направленного по нормали к поверхности S внутрь жидкости.

Уравнения (3) и (4) с граничными условиями (6) и (7) могут быть решены методом конечных элементов. Соответствующие процедуры описаны и протестированы в статьях [16, 18, 19] и ряде других работ (см. [4] и имеющиеся там ссылки). Однако, как указано во Введении, практическая реализация данного подхода требует большого объема вычислений. В следующем разделе показано, что при расчете низкочастотных компонент звукового поля объем вычислений можно существенно уменьшить. Это достигается благодаря тому, что искомое звуковое поле удается приближенно выразить через давление ${{\tilde {p}}_{f}}$ на поверхности тела S, найденное на этапе 1.

Искомое решение уравнения (3) представим в виде суперпозиции

(8)
$\tilde {p}\left( {\mathbf{x}} \right) = {{\tilde {p}}_{I}}\left( {\mathbf{x}} \right) + {{\tilde {p}}_{S}}\left( {\mathbf{x}} \right),$
где ${{\tilde {p}}_{I}}\left( {\mathbf{x}} \right)$ и ${{\tilde {p}}_{S}}\left( {\mathbf{x}} \right)$ интерпретируются как падающая на поверхность S и рассеянная на этой поверхности компоненты полного поля соответственно. Компонента ${{\tilde {p}}_{I}}\left( {\mathbf{x}} \right)$ представляет поле, которое возбудили бы источники в правой части (3) в отсутствие упругого тела. Компонента pS(x) представляет собой решение (3) с нулевой правой частью и граничными условиями на S, обеспечивающими выполнение условий (6) и (7) для суммарного поля (8).

ПАДАЮЩЕЕ ПОЛЕ

Падающее поле выражается соотношением

(9)
${{\tilde {p}}_{I}}\left( {\mathbf{x}} \right) = - \int {G\left( {{\mathbf{x}},{\mathbf{y}}} \right)} \frac{{{{\partial }^{2}}{{{\tilde {T}}}_{{ij}}}\left( {\mathbf{y}} \right)}}{{\partial {{x}_{i}}\partial {{x}_{j}}}}dV\left( {\mathbf{y}} \right),$
где dV(y) – элемент объема с центром в точке y, G(x, y) – функция Грина уравнения Гельмгольца. Поскольку источники являются акустически компактными (растояние между источниками пренебрежимо мало), функцию Грина можно записать в виде:

(10)
$G\left( {{\mathbf{x}},{\mathbf{y}}} \right) = - \frac{1}{{4{{\pi }}\left| {{\mathbf{x}} - {\mathbf{y}}} \right|}}.$

Покажем, что в этом случае падающее поле ${{\tilde {p}}_{I}}\left( {\mathbf{x}} \right)$ относительно просто выражется через давление ${{\tilde {p}}_{f}}\left( {\mathbf{x}} \right)$, найденное в приближении несжимаемой жидкости на этапе 1. Воспользуемся уравнением Эйлера несжимаемой жидкости [20]:

${{{{\rho }}}_{0}}\frac{{\partial {{v}_{i}}}}{{\partial t}} + {{{{\rho }}}_{0}}{{v}_{j}}\frac{{\partial {{v}_{i}}}}{{\partial {{x}_{j}}}} = - \frac{{\partial {{p}_{f}}}}{{\partial {{x}_{i}}}}.$

Взяв дивергенцию от обеих частей этого уравнения, т.е. вычислив производные ∂/∂xi, с учетом условия несжимаемости ∂vi/∂xi = 0 получаем

${{{{\rho }}}_{0}}\frac{{\partial {{v}_{j}}}}{{\partial {{x}_{i}}}}\frac{{\partial {{v}_{i}}}}{{\partial {{x}_{j}}}} = {{{{\rho }}}_{0}}\frac{{{{\partial }^{2}}\left( {{{v}_{i}}{{v}_{j}}} \right)}}{{\partial {{x}_{i}}\partial {{x}_{j}}}} = - \frac{{{{\partial }^{2}}{{p}_{f}}}}{{\partial x_{i}^{2}}}.$

Отсюда следует, что

$\frac{{{{\partial }^{2}}{{{\tilde {T}}}_{{ij}}}}}{{\partial {{x}_{i}}\partial {{x}_{j}}}} = - {{\Delta }}{{\tilde {p}}_{f}}.$

Подставляя это выражение в (9), находим

(11)
${{\tilde {p}}_{I}}\left( {\mathbf{x}} \right) = \int {G\left( {{\mathbf{x}},{\mathbf{y}}} \right)} {{\Delta }}{{\tilde {p}}_{f}}\left( {\mathbf{y}} \right)dV\left( {\mathbf{y}} \right).$

Воcпользуемся формулой Грина

$\begin{gathered} \int\limits_V {\left[ {{{{\tilde {p}}}_{f}}\left( {\mathbf{y}} \right){{\Delta }}G\left( {{\mathbf{x}},{\mathbf{y}}} \right) - G\left( {{\mathbf{x}},{\mathbf{y}}} \right){{\Delta }}{{{\tilde {p}}}_{f}}\left( {\mathbf{y}} \right)} \right]} dV\left( {\mathbf{y}} \right) = \\ = \int\limits_S {\left[ { - {{{\tilde {p}}}_{f}}\left( {\mathbf{y}} \right)\frac{{\partial G\left( {{\mathbf{x}},{\mathbf{y}}} \right)}}{{\partial {{y}_{i}}}} + G\left( {{\mathbf{x}},{\mathbf{y}}} \right)\frac{{\partial {{{\tilde {p}}}_{f}}\left( {\mathbf{y}} \right)}}{{\partial {{y}_{i}}}}} \right]} {\kern 1pt} {{n}_{i}}\left( {\mathbf{y}} \right)dS\left( {\mathbf{y}} \right), \\ \end{gathered} $
где V – весь объем, занятый жидкостью, dS(y) – элемент поверхности S с центром в точке y. Учитывая (11) и принимая во внимание, что внутри жидкости, а нормальная производная ${{\tilde {p}}_{f}}$ на границе S обращается в нуль, получаем искомое выражение для ${{\tilde {p}}_{I}}$ через ${{\tilde {p}}_{f}}$:

(12)
$\begin{gathered} {{{\tilde {p}}}_{I}}\left( {\mathbf{x}} \right) = {{{\tilde {p}}}_{f}}\left( {\mathbf{x}} \right) - \\ - \,\,\frac{1}{{4\pi }}\int\limits_S {{{{\tilde {p}}}_{f}}} \left( {\mathbf{y}} \right){{n}_{i}}\left( {\mathbf{y}} \right)\frac{\partial }{{\partial {{y}_{i}}}}\frac{1}{{\left| {{\mathbf{x}} - {\mathbf{y}}} \right|}}dS\left( {\mathbf{y}} \right){\kern 1pt} .{\kern 1pt} \\ \end{gathered} $

Далее при формулировке граничного условия для рассеянной компоненты ${{\tilde {p}}_{S}}$ нам понадобится нормальная производная падающей компоненты ${{\tilde {p}}_{I}}$ на границе S. Дифференцируя (12) и вновь пользуясь обращением в нуль нормальной производной ${{\tilde {p}}_{f}}$ на S, получаем выражение

(13)
$\begin{gathered} \frac{{\partial {{p}_{I}}\left( {\mathbf{x}} \right)}}{{\partial {{x}_{i}}}}{{n}_{i}}\left( {\mathbf{x}} \right) = - \frac{1}{{4\pi }}\int\limits_S {{{{\tilde {p}}}_{f}}\left( {\mathbf{y}} \right)} {\kern 1pt} {{n}_{i}}\left( {\mathbf{x}} \right){{n}_{j}}\left( {\mathbf{y}} \right) \times \\ \times \,\,\left( {\frac{{{{\partial }^{2}}}}{{\partial {{x}_{i}}\partial {{y}_{j}}}}\frac{1}{{\left| {{\mathbf{x}} - {\mathbf{y}}} \right|}}} \right)dS\left( {\mathbf{y}} \right) = \\ = \frac{{n\left( x \right)}}{{4\pi }}\int\limits_S {{{{\tilde {p}}}_{f}}} \left( {\mathbf{y}} \right)\left( {{\mathbf{n}}\left( {\mathbf{y}} \right)\nabla } \right)\nabla \frac{1}{R}dS\left( {\mathbf{y}} \right), \\ \end{gathered} $
где R = |x – y|, а оператор $\nabla {\text{\;}} = {\partial \mathord{\left/ {\vphantom {\partial {\partial {\mathbf{x}}}}} \right. \kern-0em} {\partial {\mathbf{x}}}}$, т.е. он действует только на координаты x. Мы ищем значения нормальной производной (13) в точках x, находящихся вблизи границы, но все же в толще воды. Поэтому можем воспользоваться соотношением ${{\Delta }}\frac{1}{R} = 0$. С его помощью легко проверить, что
$\left( {{\mathbf{n}}\left( {\mathbf{y}} \right)\nabla } \right)\nabla \frac{1}{R} = - \nabla \times \left( {{\mathbf{n}}\left( {\mathbf{y}} \right) \times \nabla \frac{1}{R}} \right)$
и представить (14) в виде
(14)
$\frac{{\partial {{p}_{I}}\left( {\mathbf{x}} \right)}}{{\partial {{x}_{i}}}}{{n}_{i}}\left( {\mathbf{x}} \right) = \frac{1}{{4{{\pi }}}}{\mathbf{n}}\left( {\mathbf{x}} \right)\left( {\nabla \times {\mathbf{B}}\left( {\mathbf{x}} \right)} \right),$
где

(15)
${\mathbf{B}}\left( {\mathbf{x}} \right) = - \int\limits_S {{{{\tilde {p}}}_{f}}} \left( {\mathbf{y}} \right)\left( {{\mathbf{n}}\left( {\mathbf{y}} \right) \times \nabla \frac{1}{R}} \right)dS\left( {\mathbf{y}} \right).$

Таким образом, мы получили формулу, выражающую нормальную производную падающего поля pI на границе S через ${{\tilde {p}}_{f}}$.

РАССЕЯННОЕ ПОЛЕ

Перейдем к вычислению рассеянной компоненты ${{\tilde {p}}_{S}}$, которая подчиняется уравнению Гельмгольца

(16)
${{\Delta }}{{\tilde {p}}_{S}} + {{k}^{2}}{{\tilde {p}}_{S}} = 0$
со следующим из (8) и (6) граничным условием на поверхности S

(17)
${{\left. {{{n}_{q}}\left( {\mathbf{x}} \right)\left( {\frac{{\partial {{{\tilde {p}}}_{S}}\left( {\mathbf{x}} \right)}}{{\partial {{x}_{q}}}} + \frac{{\partial {{{\tilde {p}}}_{I}}\left( {\mathbf{x}} \right)}}{{\partial {{x}_{q}}}} - {{{{\omega }}}^{2}}{{{{\rho }}}_{0}}{{{\tilde {u}}}_{q}}\left( {\mathbf{x}} \right)} \right)} \right|}_{{{\mathbf{x}} \in S}}} = 0.$

Вдали от области турбулентных пульсаций поле ${{\tilde {p}}_{S}}$ удовлетворяет условию излучения [5].

Для решений уравнений (16) и (4) воспользуемся методом конечных элементов. При этом учтем, что связь между искомыми функциями ${{\tilde {p}}_{S}}$ и ${{\tilde {u}}_{i}}$ задается граничными условиями на поверхности S: условием (17) и следующим из (7) условием

(18)
${{\left. {{{{{{\tilde {\sigma }}}}}_{{ij}}}\left( {\mathbf{x}} \right){{n}_{j}}\left( {\mathbf{x}} \right)} \right|}_{{{\mathbf{x}} \in S}}} = - {{\left. {\left[ {{{{\tilde {p}}}_{I}}\left( {\mathbf{x}} \right) + {{{\tilde {p}}}_{S}}\left( {\mathbf{x}} \right)} \right]{{n}_{i}}\left( {\mathbf{x}} \right)} \right|}_{{{\mathbf{x}} \in S}}}.$

Компоненту поля давления в жидкости ${{\tilde {p}}_{S}}\left( {\mathbf{x}} \right)$ представим в виде суперпозиции функций формы ${{{{\varphi }}}_{j}}\left( {\mathbf{x}} \right)$, $j = 1, \ldots ,J$, определенных в области V:

(19)
${{\tilde {p}}_{S}}\left( {\mathbf{x}} \right) = {{P}_{j}}{{{{\varphi }}}_{j}}\left( {\mathbf{x}} \right).$

Аналогично, компоненты вектора деформации в упругой среде ${{\tilde {u}}_{i}}\left( {\mathbf{x}} \right)$ представим в виде суперпозиции функций формы ψl(x), l = 1, …, L, определенных внутри области Ve:

(20)
${{\tilde {u}}_{i}}\left( {\mathbf{x}} \right) = {{U}_{{il}}}{{\psi }_{l}}\left( {\mathbf{x}} \right).{\text{\;\;}}$

Наша задача сводится к отысканию J неизвестных коэффициентов Pj и 3 × L коэффициентов Uil.

Умножим (16) на φj(x) и проинтегрируем полученное выражение по области V. Пользуясь теоремой Грина, находим

(21)
$\int\limits_V {\left[ { - \nabla {{{{\varphi }}}_{j}}\nabla {{{\tilde {p}}}_{S}} + {{k}^{2}}{{{{\varphi }}}_{j}}{{{\tilde {p}}}_{S}}} \right]} {\kern 1pt} dV + \int\limits_S {{{{{\varphi }}}_{j}}\frac{{\partial {{{\tilde {p}}}_{S}}}}{{\partial {{x}_{q}}}}} {\kern 1pt} {{n}_{q}}dS = 0.$

Здесь и далее для краткости опускаем аргумент x в подынтегральных выражениях. Учитывая условие (17), перепишем (21) в виде

(22)
$\begin{gathered} \int\limits_V {\left[ { - \nabla {{{{\varphi }}}_{j}}\nabla {{{\tilde {p}}}_{S}} + {{k}^{2}}{{{{\varphi }}}_{j}}{{{\tilde {p}}}_{S}}} \right]} dV = \\ = \int\limits_S {{{{{\varphi }}}_{j}}\frac{{\partial {{{\tilde {p}}}_{I}}}}{{\partial {{x}_{q}}}}} {\kern 1pt} {{n}_{q}}dS - {{{{\omega }}}^{2}}{{{{\rho }}}_{0}}\int\limits_S {{{{{\varphi }}}_{j}}{{{\tilde {u}}}_{q}}{{n}_{q}}dS} {\text{}}. \\ \end{gathered} $

Первое слагаемое в правой части (22) преобразуем с помощью (14). Получаем

(23)
$\int\limits_S {{{{{\varphi }}}_{j}}\frac{{\partial {{{\tilde {p}}}_{I}}}}{{\partial {{x}_{q}}}}} {\kern 1pt} {{n}_{q}}dS = \frac{1}{{4{{\pi }}}}\mathop \smallint \limits_S {\text{}}{{{{\varphi }}}_{j}}\left( {\mathbf{x}} \right){\mathbf{n}}\left( {\mathbf{x}} \right)\left( {\nabla \times {\mathbf{B}}} \right)dS\left( {\mathbf{x}} \right).$

С использованием формулы

$\nabla \times {{{{\varphi }}}_{j}}{\mathbf{B}} = \nabla {{{{\varphi }}}_{j}} \times {\mathbf{B}} + {{{{\varphi }}}_{j}}\nabla \times {\mathbf{B}}$
и с учетом того, что вследствие замкнутости поверхности S из теоремы Стокса [21] следует равенство
$\int\limits_S {\left[ {\nabla \times \left( {{{{{\varphi }}}_{j}}{\mathbf{B}}} \right)} \right]} {\kern 1pt} d{\mathbf{S}} = 0,$
можно получить

(24)
$\begin{gathered} \frac{1}{{4\pi }}\mathop \smallint \limits_S {\text{}}{{{{\varphi }}}_{j}}\left( {\mathbf{x}} \right){\mathbf{n}}\left( {\mathbf{x}} \right)\left( {\nabla \times {\mathbf{B}}} \right)dS\left( {\mathbf{x}} \right) = \\ = - \frac{1}{{4{{\pi }}}}\int\limits_S {\left[ {\nabla {{{{\varphi }}}_{j}}\left( {\mathbf{x}} \right) \times {\mathbf{B}}\left( {\mathbf{x}} \right)} \right]} {\kern 1pt} {\mathbf{n}}\left( {\mathbf{x}} \right)dS\left( {\mathbf{x}} \right). \\ \end{gathered} $

Подставляя найденное выражение в (22) и заменяя ${{\tilde {p}}_{s}}\left( {\mathbf{x}} \right)$ и ${{\tilde {u}}_{i}}\left( {\mathbf{x}} \right)$ суммами (19) и (20) соответственно, находим

(25)
${{Z}_{{jj{\kern 1pt} '}}}{{P}_{{j{\kern 1pt} '}}} + {{Y}_{{j,i,l}}}{{U}_{{i,l}}} = {{f}_{j}},$
где

${{Z}_{{jj{\kern 1pt} '}}} = \int\limits_V {\left[ { - \nabla {{{{\varphi }}}_{j}}\nabla {{{{\varphi }}}_{{j{\kern 1pt} '}}} + {{k}^{2}}{{{{\varphi }}}_{j}}{{{{\varphi }}}_{{j{\kern 1pt} '}}}} \right]} {\kern 1pt} dV,$
${{Y}_{{j,i,l}}} = {{{{\omega }}}^{2}}{{\rho }}\int\limits_S {{\text{}}{{{{\varphi }}}_{j}}{{n}_{i}}{{\psi }_{l}}dS,} $
${{f}_{j}} = - \frac{1}{{4\pi }}\int\limits_S {\left[ {\nabla {{{{\varphi }}}_{j}} \times {\mathbf{B}}} \right]} {\kern 1pt} {\mathbf{n}}dS.$

Таким образом, мы получили J алгебраических уравнений относительно неизвестных коэффициентов Pj и Ui,l.

Теперь обратимся к уравнению (4). Выражение (5) представим в виде

${{\tilde {\sigma }}_{{ij}}}\left( {\mathbf{x}} \right) = {{A}_{{iji{\kern 1pt} 'j{\kern 1pt} '}}}\frac{{\partial {{\psi }_{{l{\kern 1pt} '}}}\left( {\mathbf{x}} \right)}}{{\partial {{x}_{{j{\kern 1pt} '}}}}}{{U}_{{i{\kern 1pt} 'l{\kern 1pt} '}}}.$

Умножаем (4) на ψl(x) и интегрируем по Ve. Применив теорему Гаусса, находим

$ - {{\omega }^{2}}{{\rho }_{e}}\int\limits_{{{V}_{e}}} {{{{\tilde {u}}}_{i}}{{\psi }_{l}}dV} = - \int\limits_{{{V}_{e}}} {{{{\tilde {\sigma }}}_{{ij}}}} \frac{{\partial {{\psi }_{l}}}}{{\partial {{x}_{j}}}}dV + \int\limits_S {{{{\tilde {\sigma }}}_{{ij}}}{{n}_{j}}{{\psi }_{l}}dS} .$

Пользуясь (18), получаем

$\begin{gathered} - {{\omega }^{2}}{{\rho }_{e}}\int\limits_{{{V}_{e}}} {{{{\tilde {u}}}_{i}}{{\psi }_{l}}dV} = \\ = - \int\limits_{{{V}_{e}}} {{{\sigma }_{{ij}}}\frac{{\partial {{\psi }_{l}}}}{{\partial {{x}_{j}}}}dV} - \int\limits_S {\left( {{{{\tilde {p}}}_{I}} + {{{\tilde {p}}}_{S}}} \right)} {\kern 1pt} {{n}_{i}}{{\psi }_{l}}dS. \\ \end{gathered} $

Отсюда

$\begin{gathered} - {{\omega }^{2}}{{\rho }_{e}}{{U}_{{il{\kern 1pt} '}}}\int\limits_{{{V}_{e}}} {{{\psi }_{{l{\kern 1pt} '}}}{{\psi }_{l}}dV} + \int\limits_{{{V}_{e}}} {{{A}_{{iji{\kern 1pt} 'j{\kern 1pt} '}}}} \frac{{\partial {{\psi }_{{l{\kern 1pt} '}}}}}{{\partial {{x}_{{j{\kern 1pt} '}}}}}{{U}_{{i{\kern 1pt} 'l{\kern 1pt} '}}}\frac{{\partial {{\psi }_{l}}}}{{\partial {{x}_{j}}}}dV + \\ + \,\,{{P}_{j}}\int\limits_S {{{\varphi }_{j}}{{n}_{i}}{{\psi }_{l}}dS} = - \int\limits_S {{{{\tilde {p}}}_{I}}{{n}_{i}}{{\psi }_{l}}dS} . \\ \end{gathered} $

В дополнение к (25) это выражение дает еще 3L линейных алгебраических уравнений, связывающих коэффициенты Pj и Qil:

(26)
${{X}_{{ili{\kern 1pt} 'l{\kern 1pt} '}}}{{U}_{{i{\kern 1pt} 'l{\kern 1pt} '}}} + {{Q}_{{ilj}}}{{P}_{j}} = {{q}_{{il}}},$
где

${{X}_{{ili{\kern 1pt} 'l{\kern 1pt} '}}} = \int\limits_{{{V}_{e}}} {\left( { - {{\omega }^{2}}{{\rho }_{e}}{{\delta }_{{ii{\kern 1pt} '}}}{{\psi }_{{l{\kern 1pt} '}}}{{\psi }_{l}} + {{A}_{{iji{\kern 1pt} 'j{\kern 1pt} '}}}\frac{{\partial {{\psi }_{{l{\kern 1pt} '}}}}}{{\partial {{x}_{{j{\kern 1pt} '}}}}}\frac{{\partial {{\psi }_{l}}}}{{\partial {{x}_{j}}}}} \right){\kern 1pt} } dV,$
${{Q}_{{ilj}}} = \int\limits_S {{{\varphi }_{j}}{{n}_{i}}{{\psi }_{l}}dS} ,$
${{q}_{{il}}} = - \int\limits_S {{{{\tilde {p}}}_{I}}{{n}_{i}}{{\psi }_{l}}} {\kern 1pt} dS.$

В итоге для расчета неизвестных амплитуд функций формы, фигурирующих в суммах (19) и (20), мы получили систему алгебраических уравнений, заданную соотношениями (25) и (26). Решение этой системы дает компоненту звукового давления в воде pS(x) и компоненты вектора деформации в упругой среде ui(x). Компонента давления pI(x) вычисляется по формуле (12).

ВЕРИФИКАЦИЯ МЕТОДА

Верификация представленного метода проводится на примере задачи о шумоизлучении твердой квадратной пластины с длиной стороны L = 2 м и толщиной h = 0.01 м, обтекаемой потоком жидкости: скорость звука c = 1500 м/с, коэффициент кинематической вязкости ν = 1.006 × 10–6 м2/с (рис. 1). Благодаря определенному соотношению сил инерции и вязкости в потоке, обтекание задней кромки пластины сопровождается генерацией периодичных вихревых структур в виде дорожки Кармана на конкретной частоте. Соотношение сил инерции и вязкости, при котором возникает дорожка Кармана, характеризуется числом Рейнольдса Re = UL/ν, где U – скорость набегающего потока, L – характерная длина тела, ν – коэффициент кинематической вязкости. Скорости набегающего потока V1 = 1.4 м/с и V2 = 35 м/с соответствуют диапазону чисел Рейнольдса, в котором возникает дорожка Кармана для цилиндра круглого сечения [22], что также справедливо для цилиндра с прямоугольным сечением [23], и в частности для тонкой пластины с выбранными размерами.

Рис. 1.

Расчетная область гидродинамической задачи об обтекании пластины.

На рис. 1 приведена расчетная область для решения гидродинамической задачи. Расчетная область представляет собой объем жидкости в виде цилиндра с радиусом основания R = L и высотой В = 4L. Внутри выделенного объема находится квадратная пластина с указанными выше значениями параметров L и h. На входной границе 1 задана скорость набегающего потока V, нормальная к плоскости основания цилиндра. На выходной границе 2 выполняется условие обращения в нуль среднего избыточного давления. Нормальная компонента скорости жидкости на поверхности цилиндра 3 и полная скорость на поверхности пластины 4 принимаются равными нулю. Тем самым выполняются условие непротекания жидкости через боковую поверхность расчетной области и условие прилипания на пластине соответственно.

В гидродинамической части задачи движение жидкости описывается системой уравнений Навье–Стокса, осредненных по Рейнольдсу c замыканием моделью турбулентности SST Ментера [24], зарекомендовавшей себя в задачах поперечного обтекания цилиндров прямоугольного сечения [25, 26]. Жидкость считается вязкой и несжимаемой, а поверхность обтекаемого тела абсолютно твердой. Расчеты выполняются при активированной схеме High Resolution по пространству, с использованием неявной схемы интегрирования по времени второго порядка для уравнений переноса массы и импульса и первого порядка для турбулентных характеристик и двумя итерациями на одном временном шаге по рекомендациям, изложенным в работе [27].

На рис. 2б для скорости потока V2 показано пространственное распределение модуля скорости жидкости в районе задней кромки пластины в фиксированный момент времени.

Рис. 2.

К постановке гидродинамической задачи об обтекании пластины.

Полученные таким образом реализации случайных полей скорости и давления использованы для расчета компонент генерируемых акустических полей на заданных частотах  f. Расчет среднеквадратичной амплитуды давления p на поверхности сферы большого радиуса R выполнен следующими тремя способами.

1. Если размер пластины L мал по сравнению с длиной акустической волны ${{\lambda }}$ = c/f, турбулентная область, включающая твердое тело, представляет эффективный акустический диполь, ось которого перпендикулярна плоскости пластины. В любой плоскости, проходящей через эту ось, зависимость давления от угла $\theta $ между направлением на точку наблюдения и плоскостью пластины выражается соотношением

${{p}_{{\text{d}}}} = F\frac{{k\left| {\sin \theta } \right|}}{{4\pi R}},$
где F – среднеквадратичная сила, действующая со стороны турбулентных пульсаций на пластину на частоте f [5, 11]. Это выражение дает простейшую оценку возбуждаемого поля p.

2. Второй способ расчета акустического давления p базируется на прямом суммировании вкладов квадрупольных источников в правой части уравнения Лайтхилла (3). Мы использовали численную процедуру решения данной задачи, детально описанную в работе [16].

3. Наконец, в качестве третьего способа расчета давления p использован подход, развитый в данной работе.

На рис. 3 сравниваются результаты, полученные указанными тремя способами. На всех графиках представлены зависимости нормированного давления p/pd на сфере большого радиуса от угла ${{\theta }}$. Здесь показаны результаты расчета на частоте f1 = 24 Гц для скорости потока V1 (рис. 3а) и на частоте f2 = 600 Гц для скорости потока V2 (рис. 3б).

Рис. 3.

Угловая зависимость нормированного звукового давления p/pd в плоскости, проходящей через перпендикуляр к центру пластины. (а) − Скорость потока 1.4 м/с, частота 24 Гц. (б) − Скорость потока 35 м/с, частота 600 Гц.

Пунктирными линиями показаны диаграммы направленности диполя, интенсивность которого определяется вертикальной интегральной гидродинамической силой на пластине. При расчетах мы не учитываем рефракцию генерируемого турбулентностью звукового поля на течении жидкости. В монографии [5] показано, что это можно сделать путем перехода от (1) к “конвективному” уравнению Лайтхилла. Однако в наших примерах элементарные оценки показывают, что вследствие малости чисел Маха влияние рефракции на диаграммы направленности пренебрежимо мало. Сплошные черная и серая кривые на графиках получены с помощью прямого суммирования вкладов квадрупольных источников [16] и метода, предложенного в данной работе, соответственно. Расхождение между этими линиями не превышает 2% на частоте 24 Гц и 5% на частоте 600 Гц. Увеличение погрешности с частотой, по-видимому, отражает то обстоятельство, что с ростом частоты растет и отклонение точной функции Грина от используемого в нашем подходе приближенного выражения (10). На частоте 600 Гц угловая зависимость диаграммы направленности излучения существенно отличается от дипольной. Максимум излучения при этом отвечает углу 125° (направление 0° соответствует положению задней кромки пластины). Данный результат хорошо согласуется с теоретическими расчетами, представленными в работах [28, 29].

ЗАКЛЮЧЕНИЕ

Традиционный метод расчета генерации шума турбулентностью с помощью конечно-элементого моделирования [15, 16] с использованием в качестве источников тензора Лайтхилла, как правило, требует хранения и обработки большого объема данных, включающих значения поля скорости во всех узлах гидродинамической модели для всех временных отсчетов. В данной статье показано, что объем вычислений можно существенно сократить при анализе возбуждения низкочастотного звука с длиной волны, превышающей размер области турбулентных пульсаций, дающих основной вклад. Ключевая идея нашего подхода заключается в том, что при расчете компоненты звукового давления ${{\tilde {p}}_{I}}$, описывающей падающее поле, в низкочастотном приближении используется упрощенное выражение для функции Грина (10). В конечном итоге это позволяет вычислить компоненты звукового давления ${{\tilde {p}}_{I}}$ и ${{\tilde {p}}_{S}}$, а также компоненты вектора деформации ${{\tilde {u}}_{i}}$ по найденным в приближении несжимаемой жидкости значениям турбулентных пульсаций давления ${{\tilde {p}}_{f}}$ на поверхности обтекаемого тела S. Информация о пульсациях скорости при этом не используется.

Ранее близкий подход был предложен в работе [18] для описания генерации звука при обтекании абсолютно твердого тела. Он также базируется на использовании приближенного выражения (10) для функции Грина звукового поля. Наши основные результаты, выраженные соотношениями (12), (25) и (26), фактически обобщают этот подход на случай обтекания упругого тела. При построении конечно-элементной схемы расчета функций ${{\tilde {p}}_{S}}\left( x \right)$ и ${{\tilde {u}}_{i}}\left( x \right)$ важную роль сыграл вывод соотношения (14), позволившего избежать численного дифференцирования при вычислении фигурирующей в граничном условии (17) нормальной производной компоненты ${{\tilde {p}}_{I}}$.

Работа выполнена в рамках государственного задания ИПФ РАН (проект 0030-2022-0003).

Список литературы

  1. Вишняков А.Н., Макашов С.Ю. Метод оптимальной аналитической аппроксимации временных выборок в применении к анализу нестационарных периодических сигналов // Акуст. журн. 2023. Т. 69. № 2. С. 155–164.

  2. Зайцев М.Ю., Копьев В.Ф., Величко С.А., Беляев И.В. Локализация и ранжирование источников шума самолета в летных испытаниях и сравнение с акустическими измерениями крупномасштабной модели крыла // Акуст. журн. 2023. Т. 69. № 2. С. 165–176.

  3. Артельный В.В., Родионов А.А., Стуленков А.В. Повышение частотного разрешения при измерении вибраций вращающихся тел с помощью лазерной виброметрии с неподвижным лучом // Акуст. журн. 2023. Т. 69. № 3. С. 351–356.

  4. Schoder S., Kaltenbacher M. Hybrid aeroacoustic computations: state of art and new achievements // J. Theor. Comput. Acoust. 2019. V. 27. № 4. P. 1950020–1–1950020–33.

  5. Glegg S. and Devenport W. Aeroacoustics of Low Mach Number Flows. Academic Press, London, 2017. 529 p.

  6. Снегирев А.Ю. Высокопроизводительные вычисления в технической физике. Численное моделирование турбулентных течений. Санкт-Петербург, 2009. 143 с.

  7. Kajishima T., Taira K. Computational fluid dynamics: Incompressible turbulent flows. Springer, Cham, Switzerland, 2017. 364 p.

  8. Sagaut P., Deck S., Terracol M. Multiscale and multiresolution approaches in turbulence. LES, DES and hybrid RANS/LES methods: applications and guidelines. Imperial College Press, London, 2013.

  9. Lighthill M.J. On sound generated aerodynamically. I. General theory // Proc. R. Soc. Lond. 1952. V. 211(A). P. 564–587.

  10. Ffowcs J.E.W., Hawkings D.L. Sound generation by turbulence and surfaces in arbitrary motion // Proc. R. Soc. Lond. 1969. V. 264(A). P. 321–342.

  11. Howe M.S. Theory of vortex sound. Cambridge U.P., New York, 2003. 216 p.

  12. Ewert R., Schröder W. Acoustic perturbation equations based on flow decomposition via source filtering // J. Comput. Phys. 2003. V. 188. № 2. P. 365–398.

  13. Титарев В.А., Фараносов Г.А., Чернышев С.А., Батраков А.С. Численное моделирование влияния взаимного расположения винта и пилона на шум турбовинтового самолета // Акуст. журн. 2018. Т. 64. № 6. С. 737–751.

  14. Curle N. The influence of solid boundaries upon aerodynamic sound // Proc. R. Soc. Lond. 1955. V. 231(A). P. 505–514.

  15. Oberai A.A., Roknaldin F., Hughes T.J.R. Computational procedures for determining structural-acoustic response due to hydrodynamic sources // Comput. Methods Appl. Mech. Engrg. 2000. V. 190. № 3. P. 345–361.

  16. Суворов А.С., Коротин П.И., Соков Е.М. Метод конечно-элементного моделирования шумоизлучения, генерируемого неоднородностями тел, движущихся в турбулентном потоке жидкости // Акуст. журн. 2018. Т. 64. № 6. С. 756–757.

  17. Ландау Л.Д., Лифшиц Е.М. Теория упругости. М.: Наука, 1987. 248 с.

  18. Oberai A.A., Wang M. Computation of trailing edge noise from an incompressible flow calculation. Center for Turbulence Research: Proc. of the Summer Program 2000. 2000. P. 343–352.

  19. Oberai A.A., Roknaldin F., Hughes T.J.R. Computation of trailing-edge noise due to turbulent flow over an airfoil // AIAA J. 2002. V. 40. № 11. P. 2206–2216.

  20. Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М.: Наука, 1976.

  21. Корн Г., Корн Т. Справочник по математике. М.: Наука, 1974.

  22. Lienhard J.H. Synopsis of lift, drag, and vortex frequency data for rigid circular cylinders. Washington State University College of Engineering Research Division Bulletin 300. 1966. P. 1–32.

  23. Девнин С.И. Аэрогидромеханика плохообтекаемых конструкций: Справочник. Л.: Судостроение, 1983. 320 с.

  24. Menter F.R. Two-equation eddy-viscosity turbulence models for engineering applications // AIAA J. 1994. V. 32. № 8. P. 1598−1605.

  25. Tian X., Ong M.C., Yang J., Myrhaug D. Unsteady RANS simulations of flow around rectangular cylinders with different aspect ratios // Ocean Eng. 2013. V. 58. № 58. P. 208–216.

  26. Wang S., Cheng W., Du R., Wang Y. Unsteady RANS modeling of flow around two-dimensional rectangular cylinders with different side ratios at Reynolds number 6.85 × 105 // Mathematical Problems in Engineering. 2020. V. 2020. № 3. P. 1−13.

  27. Menter F.R., Sechner R. Best Practice: RANS Turbulence Modeling in Ansys CFD. Ansys Germany GmbH. Matyushenko A., NTS, St. Petersburg, Russia.

  28. Roger M., Moreau S. Back-scattering correction and further extensions of Amiet’s trailing-edge noise model. Pt. 1: theory // J. Sound Vib. 2005. V. 286. P. 477–506.

  29. Roger M., Moreau S. Back-scattering correction and further extensions of Amiet’s trailing-edge noise model. Pt. 2: application // J. Sound Vib. 2009. V. 323. P. 397–425.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Акустический журнал

Архивы выпусков Информация о журнале Отправить рукопись в журнал