舰船科学技术  2024, Vol. 46 Issue (7): 57-64    DOI: 10.3404/j.issn.1672-7649.2024.07.011   PDF    
一种基于冲激响应的瞬态噪声近似预报方法
张逸豪, 王斌, 范军, 王文欢     
上海交通大学 海洋智能装备与系统教育部重点实验室,上海 200240
摘要: 为了提高水下结构瞬态噪声的计算效率和准确性,提出一种基于冲激响应的瞬态噪声近似预报方法。该方法适用于已知结构表面法向振速时的瞬态辐射噪声预报,它通过时域快速算法求解与辐射问题相互易的脉冲点源散射问题获取结构表面各单元的冲激响应,再将各冲激响应与该处相应的法向振速信号作卷积,最后以单元面积为权重对卷积结果进行加权求和便可计算出最终的瞬态噪声场。利用2种典型的解析模型,通过与时域解析解的对比验证了该方法的正确性。针对镶嵌于圆形障板的管阵列结构,利用该近似预报方法对已确定法向振速边界条件下的瞬态辐射噪声进行计算,并结合该模型几何结构分析了考察场点处辐射噪声信号的时间特征及声能量空间指向性。
关键词: 瞬态噪声     冲激响应     互易性     时域快速算法     卷积    
An approximate prediction method of transient noise based on impulse response
ZHANG Yi-hao, WANG Bin, FAN Jun, WANG Wen-huan     
MOE Key Laboratory of Marine Intelligent Equipment and System of the Ministry of Education, Shanghai Jiaotong University, Shanghai 200240, China
Abstract: In order to improve the computational efficiency and accuracy of transient noise of underwater structures, an approximate prediction method for transient noise based on impulse response is proposed. This method is applicable to the prediction of transient radiated noise when the normal vibration velocity of a rigid structure surface is known. The impulse response of each element of the structure surface is obtained by solving a pulse point source scattering problem that is reciprocal to the radiation problem using a fast time domain algorithm, and then convolves each impulse response with the corresponding normal vibration velocity signal at that location. Finally, the transient noise field is calculated by weighted summation of the convolution results using the unit area as a weight. Using two typical analytical models, the correctness of the method is verified by comparing it with the time domain analytical solution. For a tube array structure embedded in a circular baffle, the approximate prediction method is used to calculate the transient radiated noise under a certain normal vibration velocity boundary condition. Considering the geometric structure of the model, the temporal characteristics of the radiated noise signal and the spatial directionality of the acoustic energy at the investigated field point are analyzed.
Key words: transient noise     impulse response     reciprocity     fast time domain algorithm     convolution    
0 引 言

水下辐射噪声的主要来源通常有螺旋桨噪声、机械噪声、水动力噪声3种[1]。其中,稳态机械噪声,即水下结构的稳态辐射噪声的预报方法已成熟。然而,近年来水下结构的瞬态辐射噪声也渐渐引起了国内外学者的重视。典型的水下结构瞬态噪声包括舰艇内部机器的启动、转向时的舵角变化、鱼雷发射出管等过程产生的时域上持续时间短暂的宽带噪声等。水下结构瞬态噪声因其独有的突发性、随机性、突变性等性质,利用过去处理稳态噪声的方法计算结果并不准确,亦无法描述瞬态辐射噪声的声场特性[2]。因此,对水下结构瞬态噪声的研究仍不够充分,需进行更为深入的探索研究。

目前,时域有限元法[3]、时域边界元法[4]、时域有限差分法[5]、时域有限元结合时域边界元法[6]、虚拟模态结合统计能量法[7]等作为水下瞬态结构噪声计算的重要方法,在求解时长较短的情况下可得到较为准确的收敛结果。然而,在求解过程中需对维数等离散自由度数的矩阵求逆,且由于收敛性要求数值计算的时间步长不能取得很大,计算效率受到很大制约。此外,由于有限元法和边界元法等数值方法适用于归一化波数较小时的低频噪声计算,而瞬态噪声的带宽较宽,归一化波数较大时的高频噪声计算无法使用上述数值方法得到较为准确的结果。

为了提高水下结构瞬态噪声的计算效率,特别是归一化波数较大的高频噪声计算准确性,本文提出一种基于时域冲激响应的瞬态噪声预报方法。它通过目标辐射声场与散射声场的互易定理,利用基于物理声学近似方法的时域散射声场的快速算法,求解时域冲激响应,进而对在结构表面各离散单元处时域冲激响应与法向振速的时域卷积的面积加权求和,以获得结构的瞬态辐射噪声。本文通过2种典型的解析模型验证了该方法的准确性,并利用该方法计算了复杂刚性目标的瞬态结构噪声,并分析了其时域声场的时间特性和能量的空间指向性。

1 基础理论 1.1 已知结构表面振速的时域辐射声场Kirchhoff公式

在理想流体中、小振幅声压假设前提下,基于时域格林函数的时域辐射声场的Kirchhoff公式[8]如下:

$ \begin{split} p(\stackrel{\rightharpoonup }{r},t)=\,& -\displaystyle {\int }_{-T}^{T}\displaystyle {\int }_{S}\left[G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )\frac{\partial p(\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )}{\partial n}- \right.\\[-4pt] & \left. p(\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )\frac{\partial G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )}{\partial n}\right]{\rm{d}}S{\rm{d}}\tau 。\end{split} $ (1)

式中:$ G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau ) $为时域格林函数;$ \stackrel{\rightharpoonup }{{r}^{\prime }} $为表面$ S $上的空间积分变量;$ \tau $$ - T $$ T $区间内的时间积分变量。该公式通过计算给定辐射体表面$ S $上的声压、声压的空间法向导数及声压对时间的一阶偏导数在该辐射体表面$ S $上的面积分及足够长的时间区间$ - T $$ T $内的时域积分,得到时刻$ t $位于空间点$ \stackrel{\rightharpoonup }{r} $处的时域辐射声压。

对于时域格林函数$ G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau ) $而言,它需要满足:

$ \left\{\begin{aligned} & {\nabla }^{2}G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )-\dfrac{1}{{c}^{2}}\dfrac{{\partial }^{2}G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )}{\partial {t}^{2}}=\\ & \qquad-\delta (\stackrel{\rightharpoonup }{r}-\stackrel{\rightharpoonup }{{r}^{\prime }})\delta (t-\tau ),\\ & {G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )|}_{t\leqslant \tau }=0,\\ & {\frac{\partial G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )}{\partial t}|}_{t\leqslant \tau }=0,\\ & {G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )|}_{\left|\stackrel{\rightharpoonup }{r}\right|\to \infty }=0。\end{aligned}\right. $ (2)

时域格林函数的选取方式有很多种。常见的自由场中的时域格林函数[9]为:

$ G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )=\frac{1}{4{\text{π}} R}\delta (t-\tau -\frac{R}{c})\text{ },\text{ }R=\left|\stackrel{\rightharpoonup }{r}-\stackrel{\rightharpoonup }{{r}^{\prime }}\right|。$ (3)

自由场中的时域格林函数表示$\tau $时刻位置位于$ \stackrel{\rightharpoonup }{{r}^{\prime }} $的点源在$t$时刻的场点$ \stackrel{\rightharpoonup }{r} $产生的脉冲。

对于已知结构体表面$ S $上的振速求解辐射声场的问题,可以令时域格林函数位于辐射体表面$ S $处的法向导数为0,即满足表面刚性边界条件:

$ {\frac{\partial G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )}{\partial n}|}_{\stackrel{\rightharpoonup }{r}\in S}=0。$ (4)

满足式(4)边界条件的时域格林函数,即为已知结构表面法向振速时,结构体辐射声场的时域冲激响应函数。将式(4)代入式(1),可得已知结构表面法向振速时,结构体产生的时域辐射声场:

$ p(\stackrel{\rightharpoonup }{r},t)=-{\displaystyle {\int }_{-T}^{T}{\displaystyle {\int }_{S}G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )\frac{\partial p(\stackrel{\rightharpoonup }{{r}^{\prime }},\tau )}{\partial n}{\rm{d}}S}{\rm{d}}\tau } 。$ (5)

基于该时域冲激响应函数,只需给结构体表面$ S $上声压的空间法向导数,就可表示该结构体在$t$时刻位于空间点$ \stackrel{\rightharpoonup }{r} $处的时域辐射声压。

1.2 基于互易定理的瞬态噪声近似预报方法

互易定理为线性声学中的一条普适法则。本研究基于该法则,将构建已知结构表面法向振速时的辐射声场与散射声场之间的声学互易关系。

对本文涉及的2种不同声系统进行互易性分析:一种是在声场中位置A处施加单极子声源$Q(\omega )$得到刚性结构表面位置B处的入射声压与散射声压之和${p_1}(\omega )$;另一种则是在结构表面位置B处足够小的面积微元${\rm{d}}S$施加法向振速${v_n}(\omega )$得到的声场中位置A处的辐射声压${p_2}(\omega )$。这2种结构声系统存在如下等效关系[10]

$ \frac{{{p_1}(\omega )}}{{Q(\omega )}} = \frac{{{p_2}(\omega )}}{{{v_n}(\omega ){\rm{d}}S}}。$ (6)

将式(6)写成乘积形式,并作逆傅里叶变换:

$ {{\mathcal{F}}}^{-1}[{p}_{2}(\omega )Q(\omega )]={{\mathcal{F}}}^{-1}[{p}_{1}(\omega ){v}_{n}(\omega ){\rm{d}}S]。$ (7)

由于2个时域信号卷积的傅里叶变换结果等于二者信号频谱的乘积,因此式(7)可写成:

$ {p_2}(t) \times Q(t) = {p_1}(t) \times {v_n}(t) \times {\rm{d}}S 。$ (8)

式(8)是用于描述已知结构法向振速时结构体的时域辐射问题与点源激励下满足刚性边界条件的结构体的时域散射问题互易关系式。其中,$ * $为卷积运算符。

$Q(\omega )$的幅值为${Q_0}$,则有:

$ Q(t) = {Q_0}\delta (t) 。$ (9)

即声场中位置A处的$Q(t)$是幅值为${Q_0}$的冲激脉冲信号。这样,便可利用散射问题中结构表面位置B处的入射声压与散射声压之和${p_1}(t)$与辐射问题中加载于结构表面位置B处的法向振速${v_n}(t)$进行时域卷积运算,并与面积微元${\rm{d}}S$相乘即可得到瞬态辐射问题中在声场位置A处产生的辐射声压${p_2}(t)$。即:

$ {p_2}(t) = \frac{1}{{{Q_0}}} \cdot {p_1}(t) \times {v_n}(t) \cdot {\rm{d}}S 。$ (10)

由此,可得出基于互易定理的已知结构法向振速时,结构体的时域辐射声压表达式:

$ p(\stackrel{\rightharpoonup }{r},t)=\frac{1}{{Q}_{0}}{\displaystyle {\int }_{S}{p}_{1}(\stackrel{\rightharpoonup }{{r}^{\prime }},t)\times {v}_{n}(\stackrel{\rightharpoonup }{{r}^{\prime }},t){\rm{d}}S} 。$ (11)

事实上,若交换式(5)的积分次序,时域积分项即为时域格林函数$ G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }}) $与结构表面声压的法向导数$ \dfrac{\partial p(\stackrel{\rightharpoonup }{{r}^{\prime }},t)}{\partial n} $卷积。这样,式(5)又可写成:

$ p(\stackrel{\rightharpoonup }{r},t)=-{\displaystyle {\int }_{S}G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }})\times \frac{\partial p(\stackrel{\rightharpoonup }{{r}^{\prime }},t)}{\partial n}{\rm{d}}S}。$ (12)

再由声压的法向导数与法向振速关系式:

$ \frac{\partial {v}_{n}(\stackrel{\rightharpoonup }{r},t)}{\partial t}=-\frac{1}{\rho }\frac{\partial p(\stackrel{\rightharpoonup }{r},t)}{\partial n}。$ (13)

其中:$ \rho $为水的密度。将式(13)代入式(12),并利用卷积的等效特性可得:

$ p(\stackrel{\rightharpoonup }{r},t)=\rho {\displaystyle {\int }_{S}\frac{\partial G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }})}{\partial t}\times {v}_{n}(\stackrel{\rightharpoonup }{{r}^{\prime }},t){\rm{d}}S} 。$ (14)

比较式(11)与式(14),不难发现基于点源散射与面元辐射互易关系的时域格林函数$ G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }}) $满足式(15)的条件:

$ \frac{\partial G(\stackrel{\rightharpoonup }{r},t;\stackrel{\rightharpoonup }{{r}^{\prime }})}{\partial t}=\frac{1}{\rho {Q}_{0}}{p}_{1}(\stackrel{\rightharpoonup }{{r}^{\prime }},t)。$ (15)

因此,首先利用散射问题中满足刚性边界条件的结构表面入射声压与散射声压之和$ {p}_{1}(\stackrel{\rightharpoonup }{{r}^{\prime }},t) $与辐射问题中,施加在相应结构表面位置处的法向振速$ {v}_{n}(\stackrel{\rightharpoonup }{{r}^{\prime }},t) $进行时域卷积运算,再在该结构体表面$ S $上对时域卷积结果进行面积分即可得到瞬态辐射问题中,在声场位置A处产生的辐射声压$ p(\stackrel{\rightharpoonup }{r},t) $

数值计算时,通过剖分网格的方式将目标表面离散成$N$个面单元网格,继而对每个面网格单元处的入射声压与散射声压之和与该处法向振速的时域卷积结果进行面积加权求和,最终得到离散化的时域辐射声压计算式:

$ p(\stackrel{\rightharpoonup }{r},t)=\frac{1}{{Q}_{0}}{\displaystyle \sum _{j=1}^{N}{p}_{1}(\stackrel{\rightharpoonup }{{{r}^{\prime }}_{j}},t)\times {v}_{n}(\stackrel{\rightharpoonup }{{{r}^{\prime }}_{j}},t)\Delta {S}_{j}}。$ (16)

式中:$j$为面网格单元的序号;$ \Delta {S_j} $为第$j$个面网格单元面积;$ {p}_{1}(\stackrel{\rightharpoonup }{{{r}^{\prime }}_{j}},t) $$ {v}_{n}(\stackrel{\rightharpoonup }{{{r}^{\prime }}_{j}},t) $分别为第$j$个面网格单元处的入射声压与散射声压之和及法向振速。

2 典型目标瞬态辐射声场分析

本节以单个脉动球、脉动球与刚性球阵列这2种典型模型作为研究对象,开展基于互易定理时域冲激响应的瞬态噪声近似预报方法验证。其中,脉动球与刚性球阵列模型中,包含沿直线等距排列的1个脉动球和2个刚性球,脉动球位于中间。2种典型模型中的脉动球表面均施加信号形式为单频正弦信号的法向振速。

2.1 单个脉动球的瞬态辐射声场验证与分析

该脉动球的半径为1 m。施加于球体表面的法向振速信号为单频正弦信号,频率为2 kHz,信号长度为5 ms,幅值为1 m/s。考察场点定为距球心10 m处同心球面上的某一点,不妨取为(10 m,0,0)。球表面的网格划分,以法向振速信号频率对应的1/6波长的三角形网格进行划分。图1为单个脉动球及考察场点的几何示意图。

图 1 单个脉动球及考察场点几何位置 Fig. 1 Geometric position of the single pulsating sphere and the field point

首先,进行与辐射声场相互易的散射声场的计算。在距球心10 m处的(10 m,0,0)点处设置一个单极点源,激励信号为带宽包含2 kHz的脉冲函数。通过时域快速算法求得该球表面网格每一中心点处的入射声压${p_i}(t)$与散射声压${p_s}(t)$。然后,将每一网格中心点处的入射、散射声压相加得到总声压${p_t}(t)$。接着,将每一网格中心点处的总声压${p_t}(t)$与法向振速信号${v_n}(t)$作时域卷积。最后,将每一网格对应的时域卷积结果按对应的网格面积进行加权求和,得到所求的考察点(10 m,0,0)处的时域辐射声压信号。

为了验证该近似预报算法的正确性,需将该结果与时域解析解对比。该时域解析解是基于单个脉动球的频域解析解,运用频域间接法计算得到。

图2为时域近似预报方法与时域解析解的计算结果对比图。

图 2 单个脉动球时域辐射声场的计算结果对比 Fig. 2 Comparison of calculation results of time domain radiated sound field of the single pulsating sphere

可知,二者波形吻合较好。二者考察点接收到辐射声压波形的时延一致,波形的形状与幅值略有不同。时域近似预报方法的波形幅值相对解析解较小,是由于时域快速算法中冲击响应的激励波形是以瑞克脉冲为近似波形进行计算导致的,即会有一定程度的能量泄露。此外,近似预报方法的波形在首位两端与解析解的差异是由于瑞克脉冲的波形调制导致的。

为了进一步定量比较时域近似预报方法计算结果与时域解析解计算结果的优劣,表1列出了二者所花费的计算时间以及二者时域波形的包络相关系数。

表 1 单个脉动球结果的计算时间及包络相关系数 Tab.1 Calculation time and envelope correlation coefficient of the single pulsating sphere results

可以看出,时域近似预报方法的计算时间比时域解析解要少很多,该近似方法的计算效率更高。此外,时域近似预报方法和时域解析解的包络相关系数为0.9974,相关系数接近于1,说明时域近似预报方法的结果与时域解析解的结果符合得很好。

综上,在一定的误差范围内,可以验证该时域近似预报方法的准确性和合理性。

2.2 脉动球与刚性球阵列的瞬态辐射声场验证与分析

排列于同一直线的3个球体半径均为1 m。相邻球体的球心间距为3 m。其中,对位于中心位置的脉动球表面施加法向振速信号。信号形式为单频正弦信号,频率为4 kHz,信号长度为2.5 ms,幅值为1 m/s。以中心脉动球的球心为坐标原点,以3个球的球心连线作为z轴。考察场点选取距坐标原点100 m处的同心球面上的3个点,其坐标分别为$(100,0,0)$ m、$(50\sqrt 2 ,0, 50\sqrt 2 )$ m、$(0,0,100)$ m。同样地,以法向振速信号频率对应的1/6波长作为网格最大尺寸对3个球的表面进行三角形网格划分。图3为脉动球与刚性球阵列及考察场点的几何示意图。

图 3 脉动球与刚性球阵列及考察场点几何位置 Fig. 3 Geometric positions of the array consisting of one pulsating sphere and two rigid spheres and field points

划分好网格后,首先进行与辐射声场相互易的散射声场计算。根据辐射声场考察点的数量决定了其对应散射声场的计算次数。将分别在距球心100 m处的$(100,0,0)$ m、$(50\sqrt 2 ,0,50\sqrt 2 )$ m、$(0,0,100)$ m三点处施加激励信号为带宽包含4kHz的脉冲函数单极点源。运用时域散射的快速算法分别计算3种工况下,该球表面网格每一中心点处的入射声压${p_i}(t)$与散射声压${p_s}(t)$

同样地,计算完成后将3种工况下的每一网格中心点处的入射、散射声压相加得到总声压${p_t}(t)$,并将每一网格中心点处的总声压${p_t}(t)$与法向振速信号${v_n}(t)$作时域卷积。最后,将3种工况下每一网格对应的时域卷积结果按对应的网格面积进行加权求和,分别得到所求的考察点$(100,0,0)$m、$(50\sqrt 2 ,0,50\sqrt 2 )$m、$(0,0,100)$ m处的时域辐射声压信号。

基于该脉动球与刚性球阵列的频域解析解的推导公式见附录,再利用频域间接法计算得到该脉动球与刚性球阵列模型辐射问题的时域解析解。

将以上3个考察点时域近似预报方法与解析解的计算结果进行对比。图4分别为考察点$(100,0,0)$ m、$(50\sqrt 2 ,0,50\sqrt 2 )$ m、$(0,0,100)$ m的计算结果对比图。

图 4 脉动球与刚性球阵列所有考察场点的时域辐射声场计算结果对比 Fig. 4 Comparison of time domain radiated sound field calculation results for all field points of the array consisting of one pulsating sphere and two rigid spheres

图4(a)和图4(b)可知,2个考察点的波形总体吻合较好。对每一考察点而言,时域近似预报方法得到波形相应的时延均与解析解一致。然而,近似预报方法的波形形状与幅值相较解析略有不同,特别是在二者直达接收信号的末尾处。末尾处的信号是由于中上下刚性球对脉动球辐射声场的一次散射甚至多次散射影响而产生的信号延迟效应。因为近似预报方法虽然考虑了多次散射,但未考虑到更复杂的多次散射,所以其波形与解析解的实际波形在幅值与相位上有所区别。

图4(c)中,时域近似预报方法的波形与解析解波形的直达信号拟合较好,直达信号后端毗连的波包有些差异。第1个直达波包,时域近似预报方法比解析解的幅值略低;第2个毗连波包,时域近似预报方法的信号长度小于解析解,分析是由于在求解冲激响应时受时域快速算法本身的影响,无法得到实际情况下因球体遮挡而产生的衍射效应所产生的时域波形,造成二者波形不一致。

表2列出了时域近似预报方法与时域解析解所花费的计算时间以及二者在3个考察场点时域波形的包络相关系数。

表 2 脉动球与刚性球阵列结果的计算时间及3个场点的包络相关系数 Tab.2 Calculation time and envelope correlation coefficients of three field points for the results of the array consisting of one pulsating sphere and two rigid spheres

可知,时域近似预报方法的计算时间与时域解析解相比大幅降低,该近似方法的计算效率显著提高。此外,考察场点(100, 0, 0) m、$(50\sqrt 2 ,0,50\sqrt 2 )$ m、(0, 0, 100) m的时域近似预报方法和时域解析解的包络相关系数分别为0.9376、0.9841、0.9887,三者的相关系数均接近于1,说明3个场点的2种方法得到的时域波形的包络相关系数比较理想,时域近似预报方法与时域解析解在3个场点均可认为符合得较好。

综上,在一定的误差范围内,该时域近似预报方法对脉动球与刚性球阵列的时域辐射声场实现预报的准确性得到了验证。

3 复杂结构瞬态辐射特性分析

利用该近似预报方法对结构稍复杂的模型进行时域辐射声场计算。

该复杂结构及考察场点的几何示意图如图5所示。该模型由6根完全相同且相互平行的圆柱管和一个圆形障板组成。其中,圆柱管的长度为6 m、管外直径为1 m、管厚为100 mm。圆形障板的半径为4 m、厚度为100 mm。以圆形障板的前圆端面中心作为坐标原点,并令该模型的对称平面为ZOX平面,以圆柱管的轴向方向作为X轴。a号管、b号管、c号管、d号管、e号管、f号管的轴向投影圆心分别位于(0,1 m,2 m)、(0,−1 m,2 m)、(0,3 m,0 m)、(0,1 m,0 m)、(0,−1 m,0 m)、(0,−3 m,0 m)。

图 5 复杂结构及考察场点几何位置 Fig. 5 Geometric positions of the complex structure and field points

本文将距坐标原点20 m,$x \geqslant 0$所在半空间的XY平面内的水平半圆及ZX平面内的竖直半圆上以15°为间隔的25个空间点作为辐射声压的考察场点。其中,水平半圆上的考察点,以Y轴正向作为起始角度0°,角度增大的方向为顺时针方向;竖直半圆上的考察点,以Z轴正向作为起始角度0°,角度增大的方向亦为顺时针方向。

其中,对模型中a号管的内外表面均施加法向振速信号。信号形式为单频正弦信号,频率为4 kHz,信号长度为2.5 ms,幅值为1 m/s。对该模型表面利用三角形网格进行网格划分,其中圆柱管的网格尺寸定为该法向振速信号频率对应的1/6波长,圆形障板的网格最大尺寸定为1/3波长。

3.1 近似预报方法的时域计算结果

图6为利用时域近似预报方法计算得到的水平、竖直半圆间隔15°的考察点、辐射声压时域波形对比图。

图 6 复杂结构所有考察场点时域声压计算结果 Fig. 6 Time domain sound pressure calculation results for all field points of the complex structure

可知,无论是水平方向还是竖直方向的考察场点,辐射声压信号的最早到达时刻都随着考察场点的角度变化而不同,并且其波形展宽随着角度的不同而变化。这是由于施加法向振速激励的a号管及其余结构相对于各考察场点的相对空间位置不同而造成的。

辐射声压信号的最早到达时刻取决于a号管直达辐射波的最早到达时刻和a号管的初始辐射波经模型其他结构产生的一次散射波的最早到达时刻的最小值。需注意的是,只有当考察场点与a号管之间未被模型其他结构遮挡时,a号管的直达辐射波才能被考察场点接收到。

同样地,辐射声压信号的时间延展与a号管的直达辐射波最晚到达时刻、a号管的初始辐射波,经模型其他结构产生的一次、二次甚至多次散射波最晚到达时刻有关。辐射声压信号的实际脉宽为前述两类最晚到达时刻的最大值与激励信号的脉冲宽度之和。这里,在考虑a号管直达辐射波的最晚到达时刻时,亦需判断考察场点与a号管是否被其他结构所遮挡。

3.2 声场能量的空间指向性

考察该模型在a号管施加4 kHz的单频振速信号激励下产生的辐射声场的水平方向及竖直方向的声能量空间指向性。对每一组时域声压取平方在有效时间内进行时域积分,并将归一化后的相对声能量大小作为声能量空间指向性的考察指标。图7图8分别为水平、竖直半圆间隔15°考察点辐射声能量的指向性散点图。

图 7 复杂结构辐射声能量水平方向指向性 Fig. 7 Horizontal directivity of radiated sound energy from the complex structures

图 8 复杂结构辐射声能量竖直方向指向性 Fig. 8 Vertical directivity of radiated sound energy from the complex structures

图7可知,水平方向上,在0°及15°方向的辐射声能量显著大于其他角度方向的辐射声能量。并且,在15°方向达到了最大值。是因为在0°和15°方向上,考察点位于相对于a号管的正横方向附近,且未有其他管子的遮挡影响。特别地,在180°方向上,a号管完全被与其平行的b号管遮挡,b号管的衍射效应微弱,所以其幅值接近于0。

从30°方向开始,其辐射声能量明显下降,且随着角度的能量增加呈现无规律地时增时减波动。在这些角度上,考察点远离a号管的正横方向,其振动产生的直达辐射波幅值本身较弱。另外,其5根刚性管会对a号管振动产生的辐射声场的散射效应,即由于一次散射波、二次散射波甚至多次散射波与考察点之间声程的不同,分别对考察点的直达辐射波形在时间上以不同的幅值与时延进行叠加调制。

其中,在90°方向上,辐射声能量较其他角度较大,主要是由于圆形障板对a号管振动产生的直接辐射波及其他管子的散射波散射作用,使得该考察点的辐射声能量增强。若不考虑180°方向,在45°、60°、135°以及165°方向上,其辐射声能量相对较小,这是由于管子间的互散射效应十分明显,这些方向考察点位置处的原直达辐射波的幅值被大大削减。

图8可知,竖直方向上,在0°及15°方向的辐射声能量显著大于其他角度方向的辐射声能量,并且在15°方向达到了最大值。是因为在0°和15°方向上,考察点位于相对于a号管的正横方向附近,且未有其他管子的遮挡影响。特别地,在180°方向上,a号管完全被其正下方的d号管遮挡,且d号管与e号管的衍射效应微弱,所以其声能量接近于0。

从30°方向开始,其辐射声能量明显下降,且随着角度声能量的增加呈现无规律地时增时减的波动。在这些角度上,考察点远离a号管的正横方向,其振动产生的直达辐射波幅值本身较弱。另外,其余5根刚性管及圆形障板会对a号管振动产生辐射声场的散射效应,即由于一次散射波、二次散射波甚至多次散射波与考察点之间声程的不同,分别对考察点的直达辐射波形在时间上以不同的幅值与时延进行叠加调制。

其中,在75°、90°、105°方向上,辐射声能量较其他角度较大,主要是由于圆形障板对a号管振动产生的直接辐射波及其他管子的散射波散射作用,使得该考察点的辐射声能量增强。若不考虑180°方向,在45°、135°以及165°方向上,其辐射声能量相对较小,这是由于管子与管子、管子与障板间的互散射效应十分明显,这些方向考察点位置处的原直达辐射波的幅值被大大削减。

4 结 语

1)提出一种基于时域冲激响应的瞬态噪声近似预报方法。其中,结构表面各处之于考察场点的时域冲激响应是运用时域快速算法,通过计算在考察场点设置脉冲点源激励下,结构表面各处的时域表面声压响应得到的。

2)以2种典型模型为例,通过对比该瞬态噪声近似预报方法与时域解析解的计算结果,验证了该近似预报方法的准确性,并且该近似预报方法适用于工程上对瞬态噪声的快速预报。

3)运用瞬态噪声近似预报方法,对镶嵌于圆形障板的管阵列这一复杂结构进行瞬态辐射声场的计算,并分析了所选取考察场点辐射声压信号的时间特征及声能量空间指向性。

参考文献
[1]
张勇, 张福民, 刘庆亮, 等. “深海一号”载人潜水器支持母船水下辐射噪声控制关键技术[J]. 舰船科学技术, 2022, 44(10): 49-54.
ZHANG Y, ZHANG F M, LIU Q L, et al. Research on key technologies of underwater radiated noise control of Shen Hai 1 manned submersible supporting ship[J]. Ship Science and Technology, 2022, 44(10): 49-54.
[2]
张峻铭. 水下瞬态声源声能量计算与测量方法研究[D]. 哈尔滨: 哈尔滨工程大学, 2019.
[3]
彭志. 水下目标声散射频域/时域特性仿真及应用研究[D]. 武汉: 华中科技大学, 2020.
[4]
HU F Q, PIZZO M E, NARK D M. On a time domain boudary integral equation formulation for acoustic sacttering by rigid bodies in uniform mean flow[J]. Journal of the Acoustical Socirty America, 2017, 142: 3624-3636. DOI:10.1121/1.5017734
[5]
刘佳琪, 郭永超, 朴胜春. 海洋涡旋地震成像的时域有限差分仿真[J]. 声学技术, 2022, 41(3): 426-431.
[6]
庞福振, 王洪富, 缪旭弘, 等. 瞬态冲击载荷下双层加筋圆柱壳水下声辐射研究[C]//第十八届船舶水下噪声学术讨论会论文集. 中国昆明, 2021.
PANG F Z, WANG H F, MIAO X H, et al. Study on underwater acoustic radiation of double stiffened cylindrical shells under transient impact loading[C]//Proceedings of the 18th Symposium on Ship Underwater Noise. Kunming, China, 2021.
[7]
温华兵, 吴俊杰, 马正刚, 等. 浮冰碰撞下船舶舱室瞬态噪声分析及控制[J]. 中国造船, 2020, 61(3): 121-130.
[8]
BENNETT C L, MIERAS H. Time domain integral equation solution for acoustic scattering from fluid targets[J]. Journal of the Acoustical Socirty America, 1981, 69(5): 1261-1265. DOI:10.1121/1.385808
[9]
WU S F. Transient sound radiation from impulsively accelerated bodies[J]. Journal of the Acoustical Socirty America, 1993, 94(1): 542-553. DOI:10.1121/1.407066
[10]
FAHY F J. Some applications of the reciprocity principle in experimental vibroacoustics[J]. Acoustical Physics, 2003, 49(2): 217-229. DOI:10.1134/1.1560385
[11]
GAUNAURD G C, HUANG H. Acoustic scattering by a spherical body near a plane boundary[J]. Journal of the Acoustical Socirty America, 1994, 96(4): 2526-2536. DOI:10.1121/1.410126