CN102012260B - 一种计算三维结构体流噪声声压分布的方法 - Google Patents

Abstract

本发明提出一种计算三维结构体流噪声声压分布的方法，其特征在于，首先计算三维结构体转捩区声辐射的功率谱，并将三维结构体转捩区噪声等效为分布在三维结构体表面的N个点声源；对于其中一个点声源S，计算点声源S与场点P之间沿三维结构体表面的最短距离t；计算点声源S到场点P的流噪声绕射损失，得到点声源S在场点P产生的声压为：

并积分

得到点声源S所处振动面S′在场点P处产生的声压；分别计算N个点声源所处振动面在场点P处产生的声压

根据

得到场点P处的声压G_spe。本发明能够快速方便地得到三维结构体曲面上任意点的声压大小，通过与解析解计算球面声压幅值比较，证明了本方法在整个球面上都能达到很高的计算精度。

Description

一种计算三维结构体流噪声声压分布的方法

技术领域

本发明属于水下航行体自噪声研究领域，具体为一种计算三维结构体流噪声声压分布的方法，为水下航行体等三维结构体流噪声的绕射损失计算提供方法。

背景技术

在水中航行的三维结构体表面会形成转捩区，在转捩区中会产生流噪声。目前计算流噪声在三维结构体表面的声压分布时，只有极少数问题能够求得精确解。在实际工程应用中，流噪声是通过声绕射在三维结构体表面传播，所以在计算三维结构体表面的声压分布时，大多数采用有限元或边界元等数值方法，将结构或结构表面划分为较细的网格。但有限元或边界元等数值计算方法计算规模很大，特别是针对高频计算问题，采用有限元方法计算封闭空间声模态时，由于声模态密集重叠，网格难以有效区分各个模态，导致有限元方法失效；对于自由空间的高频问题，由于边界元法的特征频率点数量大，难以完全消除各个特征频率点的影响，所以边界元方法也难以实现。

发明内容

要解决的技术问题

为了计算流噪声在三维结构体表面的声压分布，本发明提出一种计算三维结构体流噪声声压分布的方法，主要根据几何绕射理论计算转捩区产生的流噪声经过绕射后的绕射损失，从而得到三维结构体曲面上任意点的声压大小。

技术方案

所述的一种计算三维结构体流噪声声压分布的方法，其特征在于：包括以下步骤：

步骤1：确定三维结构体的流体线型，根据Lauchle提出的过渡区湍流猝发的单极子声辐射理论，计算得到三维结构体转捩区声辐射的功率谱G_sta(r，f)；并将三维结构体转捩区噪声等效为分布在三维结构体表面的N个点声源；

步骤2：对于其中一个点声源S，计算点声源S与场点P之间沿三维结构体表面的最短距离t，其中场点P为三维结构体曲面上的任一点；

步骤3：计算点声源S到场点P的流噪声绕射损失，得到点声源S在场点P产生的声压为：

其中k为流噪声的声波波数，ρ为流体介质密度，c为介质中的声速，u(s)表示点声源S点的微面积元ds的振速，V(ξ)为硬边界Fock函数，D为三维结构体表面射线的发散系数，G(kt)为Green函数；通过积分

得到点声源S所处振动面S′在场点P处产生的声压；

步骤4：重复步骤2和步骤3，分别计算出所有N个点声源所处振动面在场点P处产生的声压根据对所有声压

叠加，得到场点P处的声压G_spe，其中P_ref为参考声压值；

步骤5：重复步骤2至步骤4，计算其他场点处的声压。

有益效果

采用本发明提出的计算三维结构体流噪声声压分布的方法，根据几何绕射理论计算转捩区产生的流噪声经过绕射后的绕射损失，能够快速方便地得到三维结构体曲面上任意点的声压大小，通过与解析解计算球面声压幅值比较，证明了本方法在整个球面上都能达到很高的计算精度。

附图说明

图1：本发明的方法流程图；

图2：半球形三维结构体中选取的场点P的位置示意图；

图3：半球形三维结构体转捩区声辐射的功率谱；

图4：点声源S和场点P示意图；

图5：半球形三维结构体表面射线示意图；

图6：流噪声在1～7号场点处的声压；

图7：本发明计算结果与解析解结果比较图；

具体实施方式

下面结合具体实施例说明本发明。

实施例：

本实施例以在水中速度为19.5m/s，直径为0.48m的半球形三维结构体为例，其外形和选取的场点P点位置如图2所示。

步骤1：根据Lauchle提出的过渡区湍流猝发的单极子声辐射理论，计算得到半球形三维结构体转捩区声辐射的功率谱G_sta(r，f)，算法为：

$G_{\mathrm{sta}}(r,f)=\frac{{\mathrm{\ρ}}^2{u}_0{u}_c^2\mathrm{\Δ}{\mathrm{\δ}}^{*2}{\left(k_c\mathrm{\Δx}\right)}^2}{8{\mathrm{\π}}^2{r}^2[1+{\left(k_c\mathrm{\Δx}\right)}^2{\left(\frac{t_i{u}_c}{\mathrm{\Δx}}\right)}^2}F\left(k_c\mathrm{\Δx}\right)$

其中，系数

需要用数值积分的方法计算，本实施例中采用复合8点Gauss积分法计算；ρ为流体密度，实施例中即为水的密度；u₀为来流速度，实施例中取为19.5m/s；u_c为湍斑的迁移速度，u_c＝0.8u₀；Δδ^*为湍流边界层厚度

与层流边界层厚度

之差，

实施中为0.19mm；k_c为湍斑迁移的波数，k_c＝2πf/u_c；r为声波传播的距离，实施例中取值为1m；t_i为上升时间，实施例中取(t_iU_c)/(Δx)＝0.45，Δx为38mm；

z为无量纲量；

$N^{*}=\frac{N\left(z\right)\mathrm{\Δx}}{u_c}$

N(z)为湍斑猝发频率，取为

$N\left(z\right)=2.38\frac{u_0}{\mathrm{\Δx}}\sqrt{(1-\mathrm{\γ}\left(z\right))\ln \frac{1}{1-\mathrm{\γ}\left(z\right)}}$

γ(z)为间歇因子

$\mathrm{\γ}\left(z\right)=1-e^{-{(1+3.4z)z}^2}$

得到的半球形三维结构体转捩区声辐射的功率谱G_sta计算结果如图3所示；并将半球形三维结构体转捩区噪声等效为分布在沿半球形三维结构体表面一周的60个点声源，并且每个点声源间隔为6°。

步骤2：对于其中一个点声源S，计算点声源S与某个场点P之间沿半球形三维结构体表面的最短距离t，其中场点P为图2中所示场点中的任意一点；根据图4中所示的半球形三维结构体表面点声源S与场点P的位置关系，得到

$t=R\arccos \left(\frac{{2Z}_p{Z}_S+2\sqrt{R^2-Z_S^2}\sqrt{R^2-Z_P^2}\cos \mathrm{\φ}}{{2R}^2}\right)$

φ为点声源S和场点P在xoy平面的夹角；Z_S为点声源S的z轴坐标；Z_p为点声源P的z轴坐标；R为半球形三维结构体的半径。

步骤3：将步骤2中得到的点声源S与某个场点P之间沿半球形三维结构体表面的最短距离t代入下面公式，计算点声源S到场点P的流噪声绕射损失，得到点声源S在场点P产生的声压为：

$\mathrm{dp}=\frac{-\mathrm{jk\ρcu}\left(s\right)}{4\mathrm{\π}}\left\{2\right[1-\frac{j}{\mathrm{kt}}\left]V\left(\mathrm{\ξ}\right)\right\}\mathrm{DG}\left(\mathrm{kt}\right)\mathrm{ds}$

其中，j为虚部符号，k为流噪声的声波波数，ρ为流体介质密度，实施例中即为水的密度，c为介质中的声速，u(s)表示点声源S点的微面积元ds的振速；V(ξ)为硬边界Fock函数：

$V\left(\mathrm{\ξ}\right)=\frac{{\mathrm{\ξ}}^{\frac{1}{2}}{e}^{j\frac{\mathrm{\π}}{4}}}{2\sqrt{\mathrm{\π}}}{\∫}_{\∞e^{-j2\mathrm{\π}/3}}^{\∞}\frac{W_2^{\′}\left(\mathrm{\τ}\right)}{W_2\left(\mathrm{\τ}\right)}{e}^{-\mathrm{j\ξ\τ}}\mathrm{d\τ}$

其中，W₂(τ)为Fock型Airy函数，定义为

$W_2\left(\mathrm{\τ}\right)=\frac{1}{\sqrt{\mathrm{\π}}}{\∫}_{\∞e^{-j2\mathrm{\π}/3}}^{\∞}{e}^{\mathrm{\τZ}-Z^3/3}\mathrm{dZ}$

W′₂(τ)为W₂(τ)对τ的导数。τ和Z为复平面参数；ξ为Fock参数：

$\mathrm{\ξ}={\left(\frac{k}{2}\right)}^{\frac{1}{3}}{\∫}_S^P\frac{1}{{\left[{\mathrm{\ρ}}_g\left(t\right)\right]}^{\frac{2}{3}}}\mathrm{dt}$

ρ_g(t)为S点沿半球形三维结构体物面到P点的曲率半径；

D为三维结构体表面射线的发散系数

$D=\sqrt{\frac{\mathrm{td}{\mathrm{\ψ}}_0}{\mathrm{d\η}\left(P\right)}}=\sqrt{\frac{\mathrm{td}{\mathrm{\ψ}}_0}{{\mathrm{\ρ}}^d\mathrm{d\ψ}}}$

如图5所示，dψ₀为点声源S的两条相邻表面射线的夹角，dψ为上述两条表面射线在场点P的切线的夹角，ρ^d为表面射线束横截面的主曲率半径。对于本实施例的半球形表面：

G(kt)为Green函数：

$G\left(\mathrm{kt}\right)=\frac{e^{-\mathrm{jkt}}}{t}.$

对点声源S所处振动面S′进行积分：

$\stackrel{\‾}{p}={\∫}_{S^{\′}}\mathrm{dp}={\∫}_{S^{\′}}\frac{-\mathrm{jk\ρcu}\left(s\right)}{4\mathrm{\π}}\left\{2\right[1-\frac{j}{\mathrm{kt}}\left]V\left(\mathrm{\ξ}\right)\right\}\mathrm{DG}\left(\mathrm{kt}\right)\mathrm{ds}$

得到得到点声源S所处振动面S′在场点P处产生的声压。

对于理想刚性球面上部分球冠，其上质点沿半径方向作等幅同相谐和振动，因此点声源S所处振动面S′在场点P处产生的声压可以简化为：

$\stackrel{\‾}{p}=\frac{-\mathrm{jk\ρc}}{4\mathrm{\π}}\left\{2\right[1-\frac{j}{\mathrm{kt}}\left]V\left(\mathrm{\ξ}\right)\right\}\mathrm{DG}\left(\mathrm{kt}\right).$

步骤4：重复步骤2和步骤3，分别计算出所有60个点声源所处振动面在场点P处产生的声压

根据

对所有声压

叠加，得到场点P处的声压G_spe，其中P_ref为参考声压值，实施例中取为1μPa；

步骤5：重复步骤2至步骤4，计算图2所示其他场点处的声压。

计算结果如图6所示，图中分别给出了1至7号场点的在10kHz、15kHz、20kHz、25kHz四个频率点处的计算声压值。

图7为计算结果与和解析解计算球面声压幅值的比较，从图中可以看出，本发明的计算结果在整个半球面上都能达到很高的计算精度，即使在θ＞150度的范围内，声压幅值很小且振荡强烈，其中θ为场点P与点声源S在半球形三维结构体上过球心的半圆截面内的圆心角，本发明的计算结果仍能和解析解保持准确的一致性。

Claims (1)

1.一种计算三维结构体流噪声声压分布的方法，其特征在于：所述三维结构体为半球形三维结构体，所述计算三维结构体流噪声声压分布的方法包括以下步骤：

步骤1：确定三维结构体的流体线型，根据Lauchle提出的过渡区湍流猝发的单极子声辐射理论，计算得到三维结构体转捩区声辐射的功率谱G_sta(r，f)，算法为：

其中，系数ρ为流体密度，u₀为来流速度，u_c为湍斑的迁移速度，Δδ^*为湍流边界层厚度 

与层流边界层厚度 

之差， 

k_c为湍斑迁移的波数，k_c＝2πf/u_c，r为声波传播的距离，t_i为上升时间， 

z为无量纲量，

N(z)为湍斑猝发频率，取为

γ(z)为间歇因子

并将三维结构体转捩区噪声等效为分布在三维结构体表面的N个点声源；

步骤2：对于其中一个点声源S，计算点声源S与场点P之间沿三维结构体表面的最短距离t，其中场点P为三维结构体曲面上的任一点；

步骤3：计算点声源S到场点P的流噪声绕射损失，得到点声源S在场点P产生 的声压为：其中k为流噪声的声波波数，ρ为流体介质密度，c为介质中的声速，u(s)表示点声源S点的微面积元ds的振速，V(ξ)为硬边界Fock函数：

其中，W₂(τ)为Fock型Airy函数，定义为

W′₂(τ)为W₂(τ)对τ的导数，τ和Z为复平面参数；ξ为Fock参数：

ρ_g(t)为S点沿半球形三维结构体物面到P点的曲率半径，

D为三维结构体表面射线的发散系数，G(kt)为Green函数，其中 

通过积分

得到点声源S所处振动面S′在场点P处产生的声压；

步骤4：重复步骤2和步骤3，分别计算出所有N个点声源所处振动面在场点P处产生的声压 根据

对所有声压 

叠加，得到场点P处的声压G_spe，其中P_ref为参考声压值；

步骤5：重复步骤2至步骤4，计算其他场点处的声压。 

Images