CN110399680A - 一种浅海弹性结构辐射声场计算方法 - Google Patents

Abstract

一种浅海弹性结构辐射声场计算方法，它涉及一种浅海信道环境下弹性结构声场分布的计算方法。本发明方法：S1、设置声源；S2、获得结构声源内部等效源强度；S3、得到弹性结构在需要计算场点的辐射声场。本发明采用等效源法完整描述弹性结构的声辐射特性，声场计算结果更符合实际弹性结构辐射声场的实际情况；本发明基于表面振速反演获得结构内部等效源强度后，基于射线声学理论计算辐射声场，适用于高频段，计算速度快、效果好。

Description

一种浅海弹性结构辐射声场计算方法

技术领域

本发明涉及一种浅海信道环境下弹性结构声场分布的计算方法，基于弹性结构表面振速，联合应用等效源法与射线声学，实现该弹性结构辐射声场的快速计算，属于浅海环境下声辐射计算领域。

背景技术

水下航行器的声辐射准确预报问题具有重要的工程应用价值。一方面，辐射噪声的准确预报是确定水下航行体安全活动半径的关键，以确保自身的隐蔽性；另一方面，辐射噪声的精确仿真，也是声场特性分析和声学测量研究的重要技术手段。我国周边海洋以浅海为主，是典型的浅海波导环境。在浅海波导中，由于存在上下边界面的声反射使传统的自由场分析方法受到严重限制，所以浅海信道中的声辐射预报变得非常复杂。

浅海环境下声传播计算模型通常都是针对点声源进行建模，其经典的模型包括射线模型，简正波模型，波数积分模型以及抛物方程模型。这种将结构视作点声源的研究方法忽略了结构表面的振动分布特性和声辐射的方向性差异，仅将弹性结构视作具有相同声源级且无指向性的点声源，已经难以满足科研需求。

有限元方法也是声场计算常用的数值计算方法，但是仅适用于近场区域，对于远距离大尺度范围的声场计算，则由于网格单元数过多，计算开销太大，而不再适用。

对于弹性结构在浅海信道环境下的声辐射计算问题，在采用单一方法难以满足计算需求的情况下，一些同时采用多种计算方法的联合计算模型被开发出来；

现有联合方法：

(1)联合波叠加法与格林函数

商德江提出了波叠加法与格林函数结合的计算方法，对浅海信道下大型弹性结构声辐射做了有效的预报。这种方法应用简正波格林函数计算远场声压，而简正波格林函数在高频下阶数较大，对计算阶数的控制要求也更为严苛，计算速度也受到较大影响，所以不适宜高频计算。

(2)联合有限元法及声传播理论

安步潮联合有限元与简正波理论提出了一种计算浅海信道下复杂弹性结构声辐射的方法，并分析了其声辐射特性，类似的，钱治文提出了另外一种联合有限元和抛物方程的计算方法。两种方法都是在近场有源区需要计算近场柱面声场数据，在远场无源区域采用简正波或者抛物方程法求取完整的声场解。简正波法仅适用于低频段；抛物方程法是一种步进递推方法，远距离计算会导致误差累积，另外该方法也仅能计算某一垂直切面，而无法计算完整的三维声场。

发明内容

针对上述不足，本发明提供一种高频、快速的浅海弹性结构辐射声场计算方法。

本发明的一种浅海弹性结构辐射声场计算方法，所述方法包括：

S1、在结构声源内部设置I个点声源，作为结构声源的等效源，I个等效源的位置位于与结构声源表面共形的面上；

S2、基于弹性结构表面振速和点源格林函数，采用等效源法获得结构声源内部等效源的强度；

S3、在S2获得的等效源的强度的基础上，采用射线方法计算各个等效源的在需要计算场点的声压，将各个源的声压在相应场点叠加，得到弹性结构在需要计算场点的辐射声场。

作为优选，所述S2包括：

S21、建立柱坐标系：坐标原点在海面，海底海面边界与r轴平行，z轴垂直海面向下，结构声源的结构几何中心在z轴上，坐标为z_s；

S22、根据海面海底边界条件，求解点源的辐射声场，作为该浅海信道环境下的点源格林函数G(r,z)；

S23、基于点源格林函数G(r,z)，计算各个等效源强度{q(Q_i)}与结构声源表面振速{u(P_s)}之间的关系矩阵函数[D_si]，其中

结构声源表面有S个点振速,P_s表示第s个表面点的坐标，Q_i表示第i个等效源的坐标，i＝1,2，…I，元素D_si表示第i个等效源的源强度到第s个表面点振速的函数传递关系，表示法向导数；

S24、根据S23得到的[D_si]，结合{q(Q_i)}＝[D_si]^-1{u(P_s)}，求解各个等效源强度{q(Q_i)}。

作为优选，所述点源格林函数G(r,z)：

k_zn为深度方向波数，为第二类零阶汉克尔函数，k_rn为r方向波数，n为简正波号，N为截取的简正波模式总数，j表示虚数D为水深。

作为优选，所述S3，包括：

S31、基于射线声学理论，在需要计算场点处，求解各个特征声线：

S32、利用程函方程描述需要计算场点(r,z)处各个特征声线经过的路程：

其中，C为积分常数，n表示折射率；

S33、各个等效源的声压A：

其中，ρ为声场介质密度，c为声速，α表示声线掠射角，α₀为给定初始掠射角，表示r坐标对声源处掠射角α₀的导数，W为单位立体角内的辐射声功率，W根据{q(Q_i)}获得；

S34、将公式一和公式二带入公式三，得到各个特征声线在场点处的声压p(r,z,t)：

t表示时间；

其中，k₀为初始位置处的声波波数，k₀＝ω/c₀，ω为角频率，c₀为初始位置声速；

S35、利用S31至S34获得所有等效源在场点处的声压，叠加后得到整个弹性结构在需要计算场点处的声压。

作为优选，利用所述S1至S3依次遍历所有场点，计算得到整个辐射声场的声压分布。

作为优选，所述S31包括：

由斯奈尔定律和几何关系求解出声线的轨迹，满足以下关系：

其中，α₀为给定初始掠射角，c(z)为声速分布函数；

利用龙格库塔法数值求解公式四，获取在需要计算场点处的声线的轨迹。

本发明的有益效果，本发明联合应用了等效源法和射线法，基于弹性结构表面振速，采用等效源法获得结构内部一系列等效源的源强度，然后采用射线方法计算各个等效源的辐射声场，然后将各个源的辐射声场叠加，得到弹性结构总的辐射声场。本发明采用等效源法完整描述弹性结构的声辐射特性，声场计算结果更符合实际弹性结构辐射声场的实际情况；本发明基于表面振速反演获得结构内部等效源强度后，基于射线声学理论计算辐射声场，适用于高频段，计算速度快、效果好。

附图说明

图1为联合等效源法和射线声学的弹性结构浅海辐射声场计算方法的原理示意图；

图2为柱坐标系示意图；

图3为等效源原理示意图；

图4射线声学原理示意图；

图5为50Hz点源声场25m深度的水平分布对比图；

图6为500Hz点源声场25m深度的水平分布对比图；

图7为50Hz偶极子声场25m深度的水平分布对比图；

图8为500Hz偶极子声场25m深度的水平分布对比图；

图9为50Hz弹性结构辐射声场25m深度的水平分布对比图；

图10为300m深信道50Hz弹性结构辐射声场二维分布对比图，其中，(a)为有限元法计算结果，(b)为等效源射线法计算结果；

图11为50m深信道50Hz弹性结构辐射声场25m深度的水平分布对比图；

图12为50m深信道500Hz弹性结构辐射声场25m深度的水平分布对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

下面结合附图和具体实施例对本发明作进一步说明，但不作为本发明的限定。

本实施方式的一种浅海弹性结构辐射声场计算方法，如图1所示，联合了等效源法和射线声学。首先在弹性结构内部设置一系列的等效源，根据等效源法的原理，弹性结构的声辐射可以等效为内部各个等效源各自辐射声场的叠加，其中，等效源强度基于信道中点源格林函数和实际结构表面振速获得，等效源辐射声场有等效源强度采用射线声学理论计算获得，本实施方式包括：

S1、在结构声源内部设置I个点声源，作为结构声源的等效源，I个等效源的位置位于与结构声源表面共形的面上；

S2、基于弹性结构表面振速和点源格林函数，采用等效源法获得结构声源内部等效源的强度；

S3、在S2获得的等效源的强度的基础上，采用射线方法计算各个等效源的在需要计算场点的声压，将各个源的声压在相应场点叠加，得到弹性结构在需要计算场点的辐射声场。

优选实施例中，本实施方式的S2包括：

S21、在浅海信道环境下，建立柱坐标系，如图2所示：坐标原点在海面，海底海面边界与r轴平行，z轴垂直海面向下，由于水平方位角θ对于点源而言是周向对称的，所以在坐标系中方位角不予以考虑。结构声源的结构几何中心在z轴上，坐标为z_s；

S22、根据海面海底边界条件，求解点源的辐射声场，作为该浅海信道环境下的点源格林函数G(r,z)；

假设声速和密度只与深度z有关(水平不变)，水深为D，则点声源柱面形式的亥姆霍兹方程为，

其中p为场点声源，k＝ωc(z)为波数，ρ(z)表示位置z处的海水密度，δ为Dirac函数。

利用变量分离法，令

p(r,z)＝R(r)Z(z) (2)

将式(2)代入式(1)，解方程得到点源声场形式解

其中，n为简正波号，a_n为系数，Z_n为深度方向简正波模式函数，k_zn为深度方向波数，为第二类零阶汉克尔函数，k_rn为r方向波数；

对于一个典型情况，海面为Dirichlet边界条件，海底为经典的Neumann边界模型，密度与深度z无关，求得深度方向的简正波模式函数的具体形式，并代入式(3)，可得到该边界条件下点声源的辐射声场的声压表达式，即格林函数

其中，N为截取的简正波模式总数。

S23、基于点源格林函数G(r,z)，计算各个等效源强度{q(Q_i)}与结构声源表面振速{u(P_s)}之间的关系矩阵函数[D_si]，其中

结构声源表面有S个点振速，P_s表示第s个表面点的坐标，Q_i表示第i个等效源的坐标，i＝1,2，L I，元素D_si表示第i个等效源的源强度到第s个表面点振速的函数传递关系，表示法向导数；

S24、根据S23得到的[D_si]，结合{q(Q_i)}＝[D_si]^-1{u(P_s)}，求解各个等效源强度{q(Q_i)}；

根据等效源法理论，对于任意形状的辐射体，场点处的辐射声压是结构内部V中源的积分

其中，P,Q为全局坐标系下的格林函数场域内场点、虚拟源坐标，ρ₀为声场介质密度，ω为角频率。q(Q)为虚源源强，G(P,Q)为该边界条件下的格林函数的两点表示法。引入线性Euler公式将方程变为表面法向振速的积分表示形式

其中，表示法向导数。该式初步建立了法向振速与虚拟源强的关系。

相关研究表明，等效源面与结构表面共型，可以提升所配置的等效源的反演精度。不妨假定虚拟声源分布在结构内部一个厚度为δ_τ的虚拟球壳面σ上，则式(6)可以写为

将虚拟球壳面分成I份，每一份的表面为σ_i，如果每个球壳冠面σ_i足够小，则求和常数δ_τ可以视为常量，于是结构表面不同点P_s处法向振速近似写为

其中，q(Q_i)表示虚拟面上第i段的源强，Q_i表示第i段等效源的坐标。

若已知结构表面S个点振速{u_m}，结构内部有I个等效源{q_n}，则可将式(8)写为矩阵方程的形式

{u(P_s)}＝[D_si]{q(Q_i)} (9)

其中，元素D_si表示第i个等效源的源强度到第s个表面点振速的函数传递关系，表示为

由式(9)可得

{q(Q_i)}＝[D_si]^-1{u(P_s)} (11)

其中，[]^-1表示矩阵求逆。需要说明的是，实际应用中，矩阵[D_si]可能存在不同程度的病态性，直接求逆可能带来较大的计算误差，根据大量学者对声学逆问题的研究结果表明，采用基于L曲线的Tikhonov正则化方法会有效地抑制扰动，是在求解逆问题前有效的处理病态矩阵的手段。

按照前面对振速u(P)表达式进行离散处理的方法，对式(5)进行离散化处理，可将P点处的声压表示为I个离散点源q(Q_i)辐射声场叠加的形式，如下所示

射线声学把声波的传播看做是一束无数条垂直于等相位面的射线的传播，每一条射线与等相位面垂直，称为声线，声线途径的距离代表波传播的路程，声线经历的时间为声波传播的时间，声线束所携带的能量即为波传播的声能量。射线声学方法是在高频条件下，研究海中声传播特性的一种常用近似方法，它描述方法简洁，结果直观、清晰，在通常条件下，射线声学的数学运算也比较简洁，通常引入以下基本假定：

(1)声线的方向就是声传播的方向，声线总是垂直波振面；

(2)声线携带能量，声场中某点上的声能是所有到达该点的声线所携带能量的叠加；

(3)声线管束中能量守恒，与管外无横向能量交换。

下面给出射线声学基本的原理公式。

声波在声场介质中传播，满足波动方程

声速c＝c(r,z)，是位置的函数，设有以下声压形式解存在

其中，t表示时间，k₀＝ω/c₀为初始位置处的声波波数，ω为角频率，c₀为初始位置声速，n(r,z)＝c₀/c(r,z)为折射率，k＝k₀n(r,z)为位置(r,z)处的声波波数，A(r,z)表示声压幅值，表示声压相位，具有长度量纲，称为程函。

将式(14)代入式(13)，解得两个如下两个方程，

当远小于k²时，式(15)和(16)可得

在射线声学中，式(17)、(18)分别称作程函方程和强度方程，是射线声学的两个基本方程。

优选实施例中，本实施方式S3包括：

S31、基于射线声学理论，在需要计算场点处，求解各个特征声线：

由斯奈尔定律和几何关系可以求解出声线的轨迹，满足以下关系

其中在给定初始掠射角α₀和声速分布函数c(z)后，利用龙格库塔法数值求解此方程可以准确的完成波束追踪，获取声线的轨迹。需要说明的是，由于声线初始掠射角不同，声线在介质中的传播轨迹也不同，因此在声场中某一位置处，可能由多条声线经过，每条声线称为该位置处的一条特征声线，在使用射线理论求解声场某一位置处的声压时，需要根据式(19)求得所有特征声线。声线上点的坐标为(r，z)，z是r的函数。因此，声线轨迹坐标可以写成(r，z(r))；

S32、利用程函方程描述需要计算场点(r,z)处各个特征声线经过的路程：

其中，C为积分常数，由初始条件决定，n表示折射率；

S33、各个等效源的声压A：

由强度方程式(19)经过推导得到单根射线声强表达式

其中，W代表单位立体角内的辐射声功率，可根据等效源强度q(Q_i)求得。

根据声学基础，点源声强表达式为

其中，ρ、c为声场介质密度和声速，联合式(21)和(22)，可得该位置处的声压幅值

α表示声线在行进到距离r处的掠射角，α₀为给定初始掠射角，

S34、将式(23)和式(19)代入式(14)即可得到该条声线在位置(r,z)处的复数声压。遍历所有由某个等效源出射的经过该场点的本征声线，可以得到该等效源在该场点位置处的声压，取全部等效源叠加后即可得到整个弹性结构在该场点位置处的声压；同样地，依次遍历所有场点，即可计算得到整个辐射声场的声压分布。

实验效果对比：

(1)点源辐射声场有限元计算结果与本方法计算结果的比较

建立坐标系，坐标原点在海面，海底海面边界与r轴平行，z轴垂直海面。信道为300m深的浅海信道，海面边界为绝对软边界，海底是绝对硬边界。在信道25m深处有一50Hz的单极点源，满足的高频条件，其单极幅值为1/4πN/m。在距离点源8m的球面上，求得球面上的法向振速，使用等效源射线法求解信道中声场并与COMSOL有限元计算结果对比。在水深25m处，距离源点不同水平距离的声场分布对比如图5：

改变水深到50m，将频率增大到500Hz，按同样的步骤使用等效源射线法求解信道中声场并与COMSOL有限元计算结果对比。在水深25m处，距离源点不同水平距离的声场分布对比如图6：

图5图6可以得出，本发明所提出的方法能够有效的计算信道中点源辐射声场。

(2)偶极子源辐射声场有限元计算结果与本文方法的比较

建立坐标系，坐标原点在海面，海底海面边界与r轴平行，z轴垂直海面。信道为300m深的浅海信道，海面边界为绝对软边界，海底是绝对硬边界。在信道25m深、20m深处各有一50Hz的单极点源构成偶极子，满足的高频条件，其单极幅值分别为1/4πN/m和-1/4πN/m。在与25m处点源距离8m的球面上，求得球面上的法向振速。使用等效源射线法求解该信道中声场并与COMSOL有限元计算结果对比。在水深25m处，距离源点不同水平距离的声场分布如图7：

改变水深到50m，将频率增大到500Hz，按同样的步骤使用等效源射线法求解信道中声场并与COMSOL有限元计算结果对比。在水深25m处，距离源点不同水平距离的声场分布对比如图8：

图7图8可以得出，与点源一样，本发明所提出的方法同样能够有效的计算信道中偶极子的辐射声场。

(3)弹性结构源辐射声场有限元计算结果与本文方法的比较

建立坐标系，坐标原点在海面，海底海面边界与r轴平行，z轴垂直海面。在300m深的浅海信道中有一钢材料弹性球壳，球心在25m深处，球壳半径8m，厚度为1m。信道上表面为绝对软边界，下表面是绝对硬边界。在球壳正下方即z＝-33m处对球壳施加一点激励大小为1N，方向沿z轴指向海面，弹性球壳受激振动辐射声波。COMSOL建立模型，结构源频率满足高频条件，大小为50Hz，满足的高频条件，提取结构表面振速使用等效源射线法求解声场，并与COMSOL有限元计算结果对比。在水深25m处，距离结构中心不同水平距离的声场分布如图9。

以1m为步长，使用等效源射线法计算每个深度处声场分布，得到信道中声场的二维分布，用声压模值的分贝表示体现，并与有限元计算结果对比如图10：

信道环境不变，结构体和激励也不做任何更改，在频率不变的条件下，减少水深为50m，使用等效源射线法求解声场，并与有限元计算结果对比。在水深25m处，距离结构中心不同水平距离的声场分布如图11：

信道环境不变，结构体和激励也不做任何更改，在水深不变的条件下，增大频率为500Hz，使用等效源射线法求解声场，并与有限元计算结果对比。在水深25m处，距离结构中心不同水平距离的声场分布如图12。

由图9～12可以得出，本发明所提出的等效源射线法对弹性结构的辐射声场计算效果良好，可以作为一种新的计算方法研究相关问题。

Claims (6)

1.一种浅海弹性结构辐射声场计算方法，所述方法适用于高频段辐射声场，其特征在于，所述方法包括：

S1、在结构声源内部设置I个点声源，作为结构声源的等效源，I个等效源的位置位于与结构声源表面共形的面上；

S2、基于弹性结构表面振速和点源格林函数，采用等效源法获得结构声源内部等效源的强度；

S3、在S2获得的等效源的强度的基础上，采用射线方法计算各个等效源的在需要计算场点的声压，将各个源的声压在相应场点叠加，得到弹性结构在需要计算场点的辐射声场。

2.根据权利要求1所述的浅海弹性结构辐射声场计算方法，其特征在于，所述S2包括：

S21、建立柱坐标系：坐标原点在海面，海底海面边界与r轴平行，z轴垂直海面向下，结构声源的结构几何中心在z轴上，坐标为z_s；

S22、根据海面海底边界条件，求解点源的辐射声场，作为该浅海信道环境下的点源格林函数G(r,z)；

S23、基于点源格林函数G(r,z)，计算各个等效源强度{q(Q_i)}与结构声源表面振速{u(P_s)}之间的关系矩阵函数[D_si]，其中

结构声源表面有S个点振速,P_s表示第s个表面点的坐标，Q_i表示第i个等效源的坐标，i＝1,2，…I，元素D_si表示第i个等效源的源强度到第s个表面点振速的函数传递关系，表示法向导数；

S24、根据S23得到的[D_si]，结合{q(Q_i)}＝[D_si]^-1{u(P_s)}，求解各个等效源强度{q(Q_i)}。

3.根据权利要求1所述的浅海弹性结构辐射声场计算方法，其特征在于，所述点源格林函数G(r,z)：

k_zn为深度方向波数，为第二类零阶汉克尔函数，k_rn为r方向波数，n为简正波号，N为截取的简正波模式总数，j表示虚数D为水深。

4.根据权利要求2或3所述的浅海弹性结构辐射声场计算方法，其特征在于，所述S3，包括：

S31、基于射线声学理论，在需要计算场点处，求解各个特征声线：

S32、利用程函方程描述需要计算场点(r,z)处各个特征声线经过的路程：

其中，C为积分常数，n表示折射率；

S33、各个等效源的声压A：

其中，ρ为声场介质密度，c为声速，α表示声线掠射角，α₀为给定初始掠射角，表示r坐标对声源处掠射角α₀的导数,W为单位立体角内的辐射声功率，W根据{q(Q_i)}获得；

S34、将公式一和公式二带入公式三，得到各个特征声线在场点处的声压p(r,z,t)：

t表示时间；

其中，k₀为初始位置处的声波波数，k₀＝ω/c₀，ω为角频率，c₀为初始位置声速；

S35、利用S31至S34获得所有等效源在场点处的声压，叠加后得到整个弹性结构在需要计算场点处的声压。

5.根据权利要求4所述的浅海弹性结构辐射声场计算方法，其特征在于，利用所述S1至S3依次遍历所有场点，计算得到整个辐射声场的声压分布。

6.根据权利要求4所述的浅海弹性结构辐射声场计算方法，其特征在于，所述S31包括：

由斯奈尔定律和几何关系求解出声线的轨迹，满足以下关系：

其中，α₀为给定初始掠射角，c(z)为声速分布函数；

利用龙格库塔法数值求解公式四，获取在需要计算场点处的声线的轨迹。