Numerical Simulation of Edge Noise Using a Method Based on Synthetic Turbulence

Full Text

1. Введение

Снижение уровня подводного шума объектов морской техники является актуальной задачей. Акустическое излучение гребного винта, являющегося наиболее существенным источником шума, включает в себя тональную и широполосную составляющие. Широкополосная часть спектра шума, превалирующая в высокочастотной области, связана с взаимодействием турбулентного пограничного слоя с лопастями движителей [1]. Шум, возникающий при взаимодействии турбулентного пограничного слоя с острыми кромками тела, называется кромочным шумом. При этом большое влияние на гидродинамический кромочный шум оказывают виброакустические характеристики лопастей, что приводит к резонансному усилению излучения звука.

С развитием современных технологий появилась возможность численного прогнозирования шумов обтекания. Чаще всего используется аналогия Лайтхилла [2], концепция которой подразумевает разделение задачи на гидродинамическую и акустическую. К настоящему времени разработан ряд эффективных аналитических и численных методов решения уравнения Лайтхилла с известной правой частью, в том числе и в присутствии отражающих границ [3–5]. Аналогия Фокс Вильямса–Хокингса [5] приобрела большую популярность и широко используется при численном прогнозировании аэроакустических процессов, в том числе шумоизлучения воздушных винтов. Использован [6] метод граничных элементов с уравнениями FW–H, чтобы понять природу возникновения шума, генерируемого воздушными винтами. Определено [6] ближнее поле воздушного винта и предсказаны его акустические характеристики в дальнем поле. FW–H активно используется в работах ЦАГИ, в частности обсуждается [7, 8] методология расчета шума воздушных винтов. В ней аэродинамические характеристики винта определяются при помощи численного решения системы уравнений Эйлера на произвольных неструктурированных сетках в рамках конечно-объемного подхода, и затем для расчета акустического излучения применяется метод FW–H. Проведено [9] исследование кавитирующего и некаветирующего винтов в условиях неоднородной идеальной жидкости и показано, что прогнозы шума FW–H согласуются с результатами, полученными с помощью классического подхода с уравнениями Бернулли. Использована [10, 11] аналогия FW–H совместно с методом моделирования крупных вихрей (Large Eddy Simulation, LES) для оценки точности различных численных методологий при решении уравнения FW–H и оценки шума, создаваемого квадратным цилиндром конечного размера, погруженного в равномерный поток. Данное исследование показало, что квадрупольные источники звука являются значимыми в широком диапазоне частот. Это приводит к пониманию одного из недостатков метода FW–H и некорректности его использования в случае пренебрежения квадрупольными источниками.

Стоит отметить и несколько других недостатков аналогии FW–H. Во-первых, в ней не учитывается упругость обтекаемой конструкции, что делает аналогию неприменимой для прогнозирования гидродинамического кромочного шума. В виде исключения алгоритмы FW–H для прогнозирования подводного шума морских гребных винтов могут применяться для описания низкочастотных процессов [12, 13] в ситуациях, когда шумоизлучение вызвано дискретными лопастными частотами или кавитацией. Во-вторых, использование LES моделей при описании гидродинамических полей является дорогостоящим, и чаще пользуются гибридными численными методами. В таких методах расчетная область охватывает лишь малое пространство вблизи исследуемой неоднородности обтекаемого тела, и вихревые особенности формируются искусственно на границах перехода между зонами RANS (осредненный поток) и LES (турбулентный поток) в гибридной гидродинамической модели [14, 15]. Процессы формирования таких искусственных вихрей являются широкополосными и искажают величину и распределение нестационарных гидродинамических нагрузок, действующих на поверхность [16]. Соответственно, результаты таких расчетов могут быть некорректны.

Альтернативные FW–H подходы основаны на методе конечно-элементного моделирования шумоизлучения, например [17–19]. Акустическое излучение в этом случае рассматривается как рассеяние созданной турбулентностью псевдозвуковой волны на поверхности обтекаемого тела. Они лишены основных недостатков FW–H, однако также требуют определения флуктуаций гидродинамических полей при помощи вихреразрешающего моделирования. Моделирование турбулентных пульсаций скорости продолжает оставаться сложнейшей задачей механики сплошной среды [20, 21].

В последние годы развивается альтернативный подход, открывающий возможность гораздо менее затратного моделирования гидродинамических источников с использованием так называемой синтетической турбулентности. Термином синтетическая турбулентность (synthetic turbulence) обозначают случайные пространственно-временные распределения скорости в турбулентном течении, которые генерируются без прямого численного решения точных или приближенных уравнений гидродинамики. Входными данными для такой модели являются статистические данные о турбулентности, которые могут быть получены из менее требовательного к вычислениям CFD-моделирования, например, моделирования Навье–Стокса со средним значением Рейнольдса (Reynolds-averaged Navier–Stokes equations, RANS) [20, 21]. Алгоритмы генерации синтетической турбулентности должны удовлетворять двум ключевым требованиям: статистические характеристики получаемых распределений скорости должны быть близки к характеристикам реальной турбулентности и численная реализация алгоритма должна быть простой и не требовать значительного объема вычислений. Существуют два основных метода формирования синтетической турбулентности для анализа генерации звука, а именно метод случайных мод Фурье (RFM) и метод пространственной фильтрации.

Впервые идея такого упрощенного моделирования турбулентности, основанная на методе RFM, была предложена в работе Кречнана [22]. Метод был применен [23] для генерации стационарного турбулентного поля скоростей с использованием спектра Фон Кармана–Пао. В [24, 25] его развили, чтобы получить нестационарную турбулентную скорость в однородном конвективном среднем потоке. Затем модель RFM и линеаризованные уравнения Эйлера были объединены [26, 27] в методе генерации стохастического шума и излучения (stochastic noise generation and radiation, SNGR). Метод SNGR был улучшен [28, 29], чтобы представить конвекцию источника звука для прогнозирования шума струи.

Другая отличающаяся категория моделей [30] основана на фильтрации поля белого шума. В англоязычной литературе он известен под аббревиатурой RPM, а его ускоренный вариант обозначается FRPM. Модель основана на выборе функции корреляции турбулентных пульсаций, что неявно означает выбор соответствующего спектра и наоборот. Однако использование функции корреляции в качестве основной характеристики турбулентности позволяет гораздо проще учитывать влияние локальных неоднородностей осредненной скорости потока на генерацию звука. Поэтому при анализе генерации звука турбулентным потоком метод (F)RPM в настоящее время служит основным инструментом моделирования синтетической турбулентности. Идея этого метода реализована и развита в многочисленных работах Р. Эверта с коллегами [31–36]. Однако стоит отметить, что метод RPM больше подходит для решения задач аэроакустики. Во-первых, при больших числах Маха необходимо учитывать нелинейные эффекты, из-за чего алгоритм реализован во временной области. Во-вторых, метод также не учитывает упругость обтекаемого тела, что приемлемо для задач аэроакустики, где велико различие импедансов воздуха и обтекаемого тела.

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

2. Математическая постановка задачи

Математическое описание данного явления может выполняться с использованием уравнений Эйлера, с учетом уравнения неразрывности и уравнения состояния:

$$\begin{gathered} \frac{Du}{Dt}=-\frac{1}{{\rho }_0}\nabla p \\ \frac{\partial \rho }{\partial t}+{\rho }_0\nabla \cdot u=0 \\ \rho ={\rho }_0+\frac{p}{c^2}, \end{gathered}$$ (2.1)

где u – гидродинамическая составляющая поля скорости, p – давление, t – время, ρ – плотность жидкости, ρ₀ – условно выбранное постоянное значение плотности, c – скорость звука. При условии, что средняя скорость обтекающего тело потока существенно превышает уровни турбулентных флуктуаций скорости, а среда, в которой они движутся, акустически неоднородна (коэффициент передачи K ≠ const), источники Лайтхилла могут быть линеаризованы. Представим компоненты скорости в виде суммы средней и пульсационной составляющей: u (x, t) = U (x) + u' (x, t), где U (x) – детерминированное среднее течение, а u' (x, t) – случайная флуктуирующая компонента.

Запишем первые два уравнения системы (2.1) с учетом описанных выше замечаний:

$$\begin{gathered} \frac{\partial {u^{\text{'}}}_i}{\partial t}+U_i\nabla u^{\text{'}}+{u^{\text{'}}}_i\nabla U=-\frac{1}{{\rho }_0}\nabla p \\ \frac{\partial \rho }{\partial t}+{\rho }_0\left(\nabla \cdot U+\nabla \cdot u^{\text{'}}\right)=\mathrm{0,} \end{gathered}$$ (2.2)

где i = 1,2,3.

Предположим, что средняя скорость потока W слабо изменяется вдоль линий тока, а сами линии тока имеют малую кривизну. Введем локальную систему координат связанную с линиями тока, направленными вдоль оси z (рис. 1). Тогда U, V = 0. Средние скорости U, V, W направлены вдоль осей x, y, z соответственно. Пульсации скорости u', v', w' также направлены вдоль осей x, y, z соответственно. Поскольку рассматривается течение в пограничном слое, где градиент средней скорости потока в направлении, нормальном к поверхности обтекаемого тела, много больше, чем в направлениях вдоль и поперек потока, это позволяет пренебречь компонентами $\frac{\partial W}{\partial z}, \frac{\partial W}{\partial x}$. Применяя операцию дивергенции к первым трем уравнениям системы (2.2) и учитывая описанные замечания, можно получить следующее уравнение:

$\frac{\partial }{\partial t}\left(\nabla \cdot u^{\text{'}}\right)+2\frac{\partial W}{\partial y}\frac{\partial v^{\text{'}}}{\partial z}+W\frac{\partial }{\partial z}\left(\nabla \cdot u^{\text{'}}\right)=-\frac{{\nabla }^{2}p}{{\rho }_0}$ (2.3)

Рис. 1. Постановка задачи

С учетом второго уравнения в системе (2.2), запишем уравнение (2.3) в спектральной области, совместно с граничными условиями на рассеивающей поверхности тела и на бесконечности:

$$\begin{gathered} {\nabla }^2\tilde{p}+\frac{{\omega }^2}{c^2}\tilde{p}+\frac{i\omega }{c}M\frac{\partial \tilde{p}}{\partial z}=-2{\rho }_0\frac{\partial \widetilde{v^{\text{'}}}}{\partial z}\frac{\partial \tilde{W}}{\partial y} \\ \underset{R\to \infty }{\lim }R\left(\frac{\partial \tilde{p}}{\partial R}+ikR\right)=\mathrm{0,} \frac{\partial \tilde{p}}{\partial N}=\rho {\omega }^2{\overline{)U}}_{\Gamma }, \end{gathered}$$ (2.4)

где p – искомое давление, k = ω / c – волновое число, ω – циклическая частота, R – расстояние от источника, N – нормаль к поверхности обтекаемого тела, U – нормальное перемещение точек обтекаемого тела, Г – обтекаемая поверхность, знак «~» – Фурье-преобразование во временной области (далее этот символ опущен), i – мнимая единица, M = W/c – стремящееся к нулю число Маха, что позволяет пренебречь третьим слагаемым в первом уравнении системы (2.4).

Эту краевую задачу можно решить интегральным методом, введя K – коэффициент передачи «источник объемного ускорения – давление» (функция Грина), который характеризует излучение в заданной точке пространства из другой точки, где стоит источник, и учитывающий неоднородность пространства и упругость обтекаемого тела. Решение можно записать в интегральном виде (2.5)

$p=-2{\rho }_0\underset{V}{\int }K\frac{\partial W}{\partial y}\frac{\partial v^{\text{'}}}{\partial z}dV$, (2.5)

где V – занятая жидкостью область.

Учитывая, что в пограничном слое производная ∂W / ∂y представляет собой интенсивность осредненной величины завихренности Ω_x и для рассматриваемых случаев слабо зависит от z, в выражении (2.5) можно выполнить интегрирование по частям и получить:

$p=-2{\rho }_0\underset{V}{\int }\frac{\partial K}{\partial z}{\Omega }_x{v}^{\text{'}}dV$ (2.6)

3. Генерация случайной реализации сигнала

Анализируя выражение (2.6) можно сказать, что акустическое излучение характеризуется наличием неоднородности акустической среды, то есть неоднородность функции Грина, неоднородностью потока, то есть вязкостью и наличием пограничного слоя, и турбулентными пульсациями нормальной компоненты скорости. При этом определение функции Грина (коэффициента передачи) возможно алгоритмами [17–19] с использованием метода взаимности, а осредненную завихренность потока можно получить методами вычислительной гидродинамики с использованием недорогостоящих RANS моделей турбулентности. В выражении (2.6) вычислительную сложность представляет лишь получение реализации пульсации нормальной компоненты скорости. Предлагается ввести модель на основе синтетической турбулентности, которая генерирует реализацию пульсации нормальной компоненты скорости. В работах [31, 32] описана концепция генерации реализаций однородной случайно функции ψ (x) скалярного аргумента x путем преобразования реализации белого шума h (x) с нулевым средним ⟨h (x)⟩ = 0 и функцией корреляции

$<h\left(x^{\text{'}}\right)h\left(x^{\text{'}\text{'}}\right)>=\delta \left(x^{\text{'}}-x^{\text{'}\text{'}}\right)$, (3.1)

где угловые скобки ⟨⟩ обозначают усреднение по ансамблю реализаций.

Принимая гипотезу о замороженной турбулентности и используя соотношения из [31, 32], статистически стационарный процесс v' (x, y, z – Wt) можно представить совокупностью бегущих вдоль трубок тока гидродинамических волн со случайной комплексной амплитудой в виде:

$v^{\text{'}}=C{e}^{i\kappa z}\int {\int }_{-\infty }^{\infty }G\left(x-\epsilon \right)G\left(y-\eta \right)h\left(\epsilon , \eta \right)d\epsilon d\eta$ (3.2)

В качестве фильтра по аналогии с [28, 29] используется Гауссов фильтр:

$G\left(x, y\right)=\frac{1}{\sqrt{r_x{r}_y}}\exp \left(-\frac{\pi x^2}{2r_x^2}\right)\exp \left(-\frac{\pi y^2}{2r_y^2}\right)$, (3.3)

где r_y и r_x – радиусы корреляции пульсаций скорости v' в направлении осей y (направление градиента W) и x, и (3.2) примет вид:

$v^{\text{'}}=C{e}^{i\kappa z}\frac{1}{\sqrt{r_x{r}_y}}\int {\int }_{-\infty }^{\infty }\exp \left(-\frac{\pi {\left(x-\epsilon \right)}^2}{2r_x^2}\right)\exp \left(-\frac{\pi {\left(y-\eta \right)}^2}{2r_y^2}\right)h\left(\epsilon , \eta \right)d\epsilon d\eta$, (3.4)

где C – случайная комплексная амплитуда, κ = ω / W – гидродинамическое волновое число.

Для определения амплитуды пульсации скорости v' вычисляется спектральная плотность мощности пульсации скорости v':

$$\begin{gathered} <v^{\text{'}}{v^{\text{'}}}^{*}>=C^2\frac{1}{r_x{r}_y}\int \underset{-\infty }{\overset{\infty }{\int }}\int \underset{-\infty }{\overset{\infty }{\int }}\exp \left(-\frac{\pi {\left(x-{\epsilon }_1\right)}^2+\pi {\left(x-{\epsilon }_2\right)}^2}{2r_x^2}\right)\times \\ \times \exp \left(-\frac{\pi {\left(y-{\eta }_1\right)}^2+\pi {\left(y-{\eta }_2\right)}^2}{2r_y^2}\right)<h\left({\epsilon }_1, {\eta }_1\right)h^{*}\left({\epsilon }_2, {\eta }_2\right)>d{\epsilon }_{1}d{\eta }_{1}d{\epsilon }_{2}d{\eta }_2 \end{gathered}$$ (3.5)

Используя функцию корреляции (3.1) выражение (3.5) примет вид:

$<v^{\text{'}}{v^{\text{'}}}^{*}>=C^2\frac{1}{r_x{r}_y}\int {\int }_{-\infty }^{\infty }\exp \left(-\frac{\pi {\left(x-\epsilon \right)}^2}{r_x^2}\right)\exp \left(-\frac{\pi {\left(y-\eta \right)}^2}{r_y^2}\right)d\epsilon d\eta$ (3.6)

Пользуясь интегралом Гаусса ${\int }_{-\infty }^{\infty }\exp (-a{x}^2)dx=-\sqrt{\pi /a}$ (a > 0) для интегрирования выражения (3.6), получим:

$C=\overline{v^{\text{'}}}$, (3.7)

где $\overline{v^{\text{'}}}=\sqrt{<v^{\text{'}}{v^{\text{'}}}^{*}>}$.

Окончательно (3.2) примет вид:

$v^{\text{'}}=\overline{v^{\text{'}}}{e}^{i\kappa z}\frac{1}{\sqrt{r_x{r}_y}}\int {\int }_{-\infty }^{\infty }\exp \left(-\frac{\pi {\left(x-\epsilon \right)}^2}{2r_x^2}\right)\exp \left(-\frac{\pi {\left(y-\eta \right)}^2}{2r_y^2}\right)h\left(\epsilon , \eta \right)d\epsilon d\eta$ (3.8)

С другой стороны, для спектральной плотности мощности пульсаций скорости можно использовать спектр турбулентности фон Кармана [37–39]:

$\overline{v^{\text{'}}}=\frac{A}{{\left(1+B\right)}^{11/12}}\sqrt{1+\frac{8}{3}B}$ (3.9)

$A=\sqrt{\frac{4E}{3\pi f_0}}$ (3.10)

$B={\left(2.678\frac{f}{f_0}\right)}^2$ (3.11)

$f_0=\frac{W}{\pi \left(r_x+r_y\right)}$, (3.12)

где E – кинетическая энергия турбулентности, f – частота. Интегрирование квадрата выражения (3.9) в частотной области приводит к кинетической энергии турбулентности, откуда может быть получено определение для константы A (3.10). Зависимость B (3.11) представлена эмпирической зависимостью [37]. Частота f₀ является эмпирической зависимостью, полученной в ходе численных экспериментов.

Таким образом (2.6) можно переписать в виде:

$p=-2{\rho }_0\underset{z}{\int }\underset{s}{\int }h\left(\epsilon , \eta \right)\underset{S}{\int }\frac{\partial K}{\partial z}\overline{v^{\text{'}}}{e}^{i\kappa z}\frac{{\Omega }_x}{\sqrt{r_x{r}_y}}\exp \left(-\frac{\pi {\left(x-\epsilon \right)}^2}{2r_x^2}\right)\exp \left(-\frac{\pi {\left(y-\eta \right)}^2}{2r_y^2}\right)dsdSdz$, (3.13)

где s – поверхность перпендикулярная линиям тока по аналогии с S.

В ходе многочисленных численных экспериментов выяснено, что радиусы корреляции r_x и r_y являются в достаточной степени малыми, что позволяет выполнить упрощенное вычисление внутреннего интеграла в (3.13) и представить звуковое давление в виде:

$p=-4{\rho }_0{\int }_z{\int }_{S}h\left(\epsilon , \eta \right){\Omega }_x\frac{\partial K}{\partial z}{e}^{i\kappa z}\overline{v^{\text{'}}}\sqrt{r_x{r}_x}dSdz$ (3.14)

Исходя из (3.14) и того, что h является дельта-коррелированной случайной функций можно записать следующее выражение для мощности звукового давления:

$\left|p\right|=4{\rho }_0\sqrt{{\int }_S{\left|{\int }_z{\Omega }_x\frac{\partial K}{\partial z}{e}^{i\kappa z}\overline{v^{\text{'}}}\sqrt{r_x{r}_y}dz\right|}^{2}dS}$ (3.15)

Анализируя физический смысл выражения (3.14) можно сказать, что звуковое поле (3.15) интерпретируется в виде некогерентной вклада в поле турбулентных пульсаций в элементарных струйках тока площадью dS. При этом существенную роль на спектральные характеристики шума оказывает пространственное Фурье – преобразование производной коэффициента передачи K по направлению линии тока.

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

Для обеспечения упрощения вычисления в (3.15) несобственного интеграла вдоль линии тока в подынтегральное выражение можно внести пространственный Гауссов фильтр F, размер которого существенно превышает область больших градиентов величины K и ограничить область интегрирования отрезком линии тока длиной l. Кроме того, внешний интеграл можно привести к дискретному виду и окончательно получить:

${\left|p\right|}^2=4{\rho }_0\sqrt{{\sum }_j{\left|{\int }_{l_j}{\Omega }_{xj}\frac{\partial K_j}{\partial z}{e}^{i{\kappa }_{j}l}\overline{{v^{\text{'}}}_j}\sqrt{r_{xj}{r}_{yj}}d{l}_j\right|}^{2}d{S}_j}$ (3.16)

Интегрирование в (3.16) выполняется по заданной сетке линий тока с заданным шагом dl вдоль линии тока. Для радиусов корреляции можно ввести модельные зависимости. Например, радиусы корреляции в приложении к задачам обтекания кромок лопастей движителя можно приблизительно взять в следующем виде:

$r_{xj}=0.2{\delta }_{0.95j}, r_{yj}=\frac{y_j}{2}$, (3.17)

где δ_0.95 – толщина пограничного слоя, определяемая из численного расчета по уровню 95% от скорости W_∞. Подробнее о получении этих модельных зависимостей изложено в разд. 5.

4. Численная реализация алгоритма

Входными данными для расчета описанным методом является RANS расчет в исследуемом объеме жидкости, включающим обтекаемое тело, и коэффициент передачи (КП) «источник объемного ускорения – давление». Определение КП проводится по методу взаимности на конечно-элементной сетке, согласно методике [17], источник объемного ускорения располагается в точке, где проводится оценка шумоизлучения. Конечно-элементная (КЭ) сетка должна включать в себя водные и структурные элементы, а также содержать набор узлов начальной линии, от которой будет начинаться построение линий тока.

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

После этого вычисляется значение интеграла по формуле (3.16). Также пользователь самостоятельно задает модельные зависимости радиусов корреляции.

5. Верификация

Верификация представленного метода проводится на примере задачи о шумоизлучении задней кромки лопасти, обтекаемой потоком жидкости: скорость звука c = 1500 м/с, коэффициент кинематической вязкости ν = 1.006 × 10^–6 м²/с. Для возможности верификации по методике, описанной в [17], рассматривается не вся кромка целиком, а несколько фрагментов, вырезанных по радиусу лопасти. То есть рассматривается набор квазидвумерных задач, отличающихся друг от друга углом установки профиля, толщиной профиля и гидродинамическими характеристиками набегающего потока.

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

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

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

На рис. 3 показано пространственное распределение модуля скорости жидкости, кинетической энергии турбулентности и модуля завихренности в районе задней кромки лопасти в фиксированный момент времени.

Рис. 3. Осредненное поле модуля завихренности (а), кинетической энергии турбулентности (б) и модуля скорости (в) в сечении вдоль потока

На рис. 4 показана КЭМ используемая для определения КП в соответствии с [42] в связной постановке для СЛАУ относительно неизвестных величин потенциала перемещений в жидкости и перемещений в упругом теле. КЭМ представляет собой объем жидкости, ограниченный сферической поверхностью, внутри которой находится исследуемая лопасть целиком. На сферической границе расчетной области специальные КЭ формируют импедансное условие не отражения проходящих по нормали через нее акустических волн, обеспечивая излучение звука по аналогии с бесконечным свободным пространством. На поверхности лопасти задано граничное условие Fluid Solid Interface (FSI), которое связывает степени свободы перемещения в структурных элементах и давления в жидкостных элементах модели. Корневое сечение лопасти жестко зафиксировано по перемещениям. Размер элемента вблизи профиля выбирается исходя из минимизации погрешности интерполяции гидродинамических полей на КЭ сетки по соотношению l_ак ≤ 4l_гд, где l_ак размер элемента в акустической КЭМ, а l_гд – размер ячейки в гидродинамической сеточной модели. В остальной области размер элемента соответствует l_ак1 ≤ λ_min / 3, где λ_min = с / f_max – наименьшая длина волны, соответствующая максимальной моделируемой частоте, c – скорость распространения звука. Используются квадратичные элементы акустической численной модели для обеспечения достаточной точности вычисления производной КП в (3.16). Количество конечных элементов АКМ составляет ~ 2 млн. Источник объемного ускорения единичной амплитуды необходимый для расчета КП методом взаимности [43] располагается в дальней зоне исследуемого объекта. На рис. 5 показана визуализация КП в сечении поперек размаха лопасти, темный цвет соответствует отрицательным значениям КП, светлый – положительным.

Рис. 4. Конечно-элементная модель для определения коэффициента передачи

Рис. 5. Реальная часть коэффициента передачи на частоте δ_0.95f / W_∞ = 0.4

На рис. 6, а показан радиус корреляции в направлении вдоль радиуса лопасти, на рис. 6,б показан радиус корреляции в нормальном направлении к поверхности лопасти на различных частотах. Серые сплошные линии соответствуют частотам $\frac{{\delta }_{0.95}f}{W_{\infty }}<0.2$, серые пунктирные – $0.2>\frac{{\delta }_{0.95}f}{W_{\infty }}>0.4$, серые штрихпунктирные – $0.4>\frac{{\delta }_{0.95}f}{W_{\infty }}>0.7$, черная – усредненная. Предложена модель, согласно которой осредненный радиус корреляции в направлении вдоль радиуса лопасти равен 20% от толщины пограничного слоя r_x = 0.2δ_0.95, а осредненный радиус корреляции в нормальном направлении к поверхности лопасти равен r_y = 2y / 3. Модельные зависимости показаны черными пунктирными линиями на рис. 4. Можно заметить, что в случае с радиусом корреляции в направлении вдоль задней кромки лопасти модельная зависимость не совпадает с усредненным радиусом корреляции вблизи поверхности лопасти и больше схожа на модельную зависимость для радиуса корреляции в нормальном направлении к поверхности лопасти. Поэтому другой моделью описания осредненного радиуса корреляции в направлении вдоль радиуса лопасти можно дополнительно рассмотреть зависимость r_x = r_y = 2y / 3.

Рис. 6. Радиус корреляции в направлении вдоль радиуса лопасти (а) и радиус корреляции в нормальном направлении к поверхности лопасти (б) на различных частотах

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

(i) Первый способ расчета акустического давления p базируется на прямом суммировании вкладов квадрупольных источников в правой части уравнения Лайтхилла. Используется численная процедура решения данной задачи, детально описанная в работе [17].

(ii) В качестве второго способа расчета давления p используется подход, развитый в данной работе.

На примере одного из фрагментов на рис. 7 показаны результаты исследования модельной зависимости осредненного радиуса корреляции в направлении вдоль радиуса лопасти. Средняя погрешность во всем рассматриваемом частотном диапазоне в случае использования модели r_x = 0.2δ составляет 3 дБ, в случае использования модели r_x = 2y / 3 – 4 дБ. Однако можно видеть, что в диапазоне частот 0.4 > δ_0.95f / W_∞ > 1.3 использование модели r_x = 2y / 3 дает заниженный результат в среднем на 5 дБ, тогда как модель r_x = 0.2δ_0.95 завышает результат в среднем на 2 дБ в диапазоне частот до δ_0.95f / W_∞ > 1.3. Таким образом, наиболее близкий результат возможно получить, используя модельную зависимость r_x = 0.2δ_0.95. Во всех представленных расчетах используется эта модель для описания осредненного радиуса корреляции в направлении вдоль радиуса лопасти.

Рис. 7. Зависимость оценки КШ от модели описания осредненного радиуса корреляции в направлении вдоль радиуса лопасти

На рис. 8–10 показана сходимость метода в зависимости от заданных пользователем параметров дискретизации линий тока – шаг вдоль начальной линии da, шаг вдоль линии тока dl, количество точек по нормали N_δ, количество точек вдоль линии тока N_z и коэффициента геометрической прогрессии к поверхности лопасти r. На рис. 8 показана сходимость метода в зависимости от длины линии тока, которая может быть представлена как N_zdl, и количества расчетных точек вдоль нее, то есть шага dl. В рамках этого исследования на толщину пограничного слоя приходится 20 слоев расчетных точек (N_δ = 20) со сгущением к поверхности с коэффициентом геометрической прогрессии r = 1.2, вдоль задней кромки 10 слоев на толщину пограничного слоя (da = 1 мм). Увеличение уровня на высоких частотах  объясняется в одном случае слишком редким шагом вдоль линии тока, а в другом – слишком маленькой длиной линии тока. Таким образом, на длину линии тока должно укладываться хотя бы 14 длин гидродинамических волн, а линия тока должна быть разбита как минимум на 50 частей.

Рис. 8. Исследование сходимости метода от параметров «количество точек вдоль линии тока» и «шаг вдоль линии тока»

Рис. 9. Исследование сходимости метода от параметра «прогрессия по нормали к поверхности»

Рис. 10. Исследование сходимости метода от параметра «количество точек по нормали к поверхности» и «шаг вдоль начальной линии»

На рис. 9 показана зависимость результата от коэффициента геометрической прогрессии. Можно видеть, что он не влияет на результат. Для построения этих зависимостей на толщину пограничного слоя приходится 20 слоев расчетных точек (N_δ = 20) со сгущением к поверхности с коэффициентом геометрической прогрессии r = 1.2, вдоль задней кромки 10 слоев на толщину пограничного слоя (da = 1 мм), а также соблюдаются требования к длине линии тока и ее дискретизации. Таким образом, коэффициент геометрической прогрессии к поверхности лопасти r можно выбрать произвольно, однако для дальнейших верификационных расчетов принято r = 1.2, поскольку этот коэффициент сгущения к стенке рекомендован при построении сеточных гидродинамических моделей для описания пограничного слоя [44].

На рис. 10 показана зависимость результата от количества точек по нормали к поверхности N_δ и шага вдоль начальной линии da. Можно видеть, что на толщину пограничного слоя достаточно двадцати слоев, а вдоль начальной линии достаточно четырех слоев на толщину пограничного слоя. Другие параметры расчета удовлетворяют требованиям, описанным выше.

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

Рис. 11. Вклад источников гидродинамического шума по нормали к поверхности лопасти

На рис. 12 сравниваются результаты, полученные указанными двумя способами. На графике представлены зависимости звукового давления в точке на сфере большого радиуса. Сплошные линии соответствуют расчетам прямым методом [17], пунктирные линии получены с помощью метода, предложенного в данной работе. Фрагмент № 1 соответствует торцевой части лопасти, фрагменты № 2 и № 3 располагаются ниже по радиусу соответственно. Для фрагмента № 1 средняя погрешность составляет 4 дБ, для фрагмента № 2 – 1 дБ, для фрагмента № 3 – 3 дБ. Можно при этом отметить качественное совпадение результатов, полученных двумя методами. Приведенные верификационные результаты, демонстрирующие высокую точность разработанного метода подтверждаются валидацией, выполненной по известным экспериментальным данным, с которыми можно ознакомиться в [45].

Рис. 12. Звуковое давление в точке на контрольной дистанции от поверхности лопасти

Заключение

В работе представлен и верифицирован метод моделирования шумоизлучения упругих тел методом на основе синтетической турбулентности. Пульсации скорости генерируются на основе информации об осредненных гидродинамических полях, путем пространственной фильтрации белого шума с заданными корреляционными характеристиками. На примере гребного винта показано, что метод позволяет прогнозировать излучение звука с низкой погрешностью, не превышающей несколько дБ. Метод реализован в ПО «САТЕС» и может использоваться в практике акустического проектирования изделий с важностью вклада гидродинамических шумов в широком частотном диапазоне.

Работа выполнена в рамках Государственного задания Минобрнауки России № FFUF-2022-0003.

About the authors

N. V. Balakireva

Author for correspondence.
Email: balakireva@ipfran.ru
Russian Federation, Nizhny Novgorod

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Problem statement

Download (11KB)

Indexing metadata

3. Fig. 2. The DG computational domain of the problem on the flow of a fragment of the blade trailing edge

Download (10KB)

Indexing metadata

4. Fig. 3. Averaged field of the curl modulus (a), turbulence kinetic energy (b) and velocity modulus (c) in a cross section along the flow path

Download (10KB)

Indexing metadata

5. Fig. 4. Finite element model for determining the transmission coefficient

Download (50KB)

Indexing metadata

6. Fig. 5. Real part of the transmission coefficient at frequency δ0.95f / W∞ = 0.4

Download (17KB)

Indexing metadata

7. Fig. 6. Correlation radius in the direction along the blade radius (a) and correlation radius normal to the blade surface (b) at different frequencies

Download (27KB)

Indexing metadata

8. Fig. 7. Dependence of the CS estimate on the model description of the averaged correlation radius in the direction along the blade radius

Download (22KB)

Indexing metadata

9. Fig. 8. Study of convergence of the method from the parameters ‘number of points along the current line’ and ‘step along the current line’

Download (29KB)

Indexing metadata

10. Fig. 9. Study of convergence of the method from the parameter ‘progression along the normal to the surface’

Download (17KB)

Indexing metadata

11. Fig. 10. Study of convergence of the method from the parameter ‘number of points along the normal to the surface’ and ‘step along the initial line’

Download (24KB)

Indexing metadata

12. Fig. 11. Contribution of hydrodynamic noise sources along the normal to the blade surface

Download (15KB)

Indexing metadata

13. Fig. 12. Sound pressure at a point at the reference distance from the blade surface

Download (29KB)

Indexing metadata