CN112257311A - Pekeris波导下结构声振计算的FE/WSM方法 - Google Patents

Abstract

本发明公开了一种Pekeris波导下结构声振计算的FE/WSM方法，包括：采用有限元法建立局域Pekeris波导下弹性结构声辐射多物理场耦合模型，进行三维Pekeris波导下弹性结构离散处理；通过有限差分法联立求解各个耦合模型，提取局部Pekeris波导下结构表面振动数据矩阵；采用波叠加法进行弹性结构内部虚拟点声源的空间优化布放，并基于球面波分解的精准镜像源法求解近场Green函数，获取各个虚拟源与结构表面离散节点的传递矩阵，通过振动数据、传递矩阵的逆运算求解近场各个虚拟源的声源强度；通过简正波扰动理论求解远场Green函数，以及计算各个虚拟源到场点的单极矩阵；将声源强度与对应的单极矩阵相乘获取声场；将各个离散虚拟源在远场场点产生的声场进行叠加求和，以等效为连续弹性结构的辐射声场。

Description

Pekeris波导下结构声振计算的FE/WSM方法

技术领域

本发明涉及浅海环境下弹性结构声辐射预报、声学测量、水中目标探测和识别等领域，尤其涉及一种基于FE/WSM(有限元/波叠加法)的Pekeris(皮克里斯)波导下弹性结构声振计算方法。

背景技术

由于人类对资源的开发和保护活动中逐渐由陆地转移到海洋，在由大陆架向海洋延伸的广泛过渡区域，属于水深小于200米的典型浅海环境，在该区域存在大量的水下航行器、输送管道和海洋平台等弹性结构，弹性结构在浅海下的声振特性研究近年吸引了学者的极大关注。对海洋波导中辐射噪声的实时监测预报、噪声控制和目标探测识别具有举足轻重的作用，是以后我国水声技术领域长期关注的热点和难点问题之一。

然而，目前弹性结构在海洋波导下的振动与声辐射研究的流体域大多考虑为无界或半空间流体域，对于海洋波导下弹性结构多物理场耦合声辐射问题的研究尚不多见。不同于以往自由场环境下弹性结构声辐射机理，浅海下的结构耦合振动和辐射声场会受到水面和水底边界反射声的重要影响，首先边界反射声会作为激励条件反作用于结构表面和附近声场，改变结构的激励属性和流固耦合作用，其次，边界反射声会与直达声干涉叠加远场声场，使得辐射声场能够在上下边界的限制下远距离传播。所以，其辐射声场计算是一个复杂的物理问题，需要充足的物理模型来描述结构与周围的流固耦合、边界反射声与结构的声壳耦合和在海底交互作用面上声与海底介质的耦合作用，而且计算方法不仅需要能分析上下边界影响下的近场耦合振动和辐射声场特性；还能研究辐射声场在浅海中的传播规律。

波导下结构声场问题涉及复杂的多物理场耦合环境，在这些耦合条件下数学理论推导难以求解，数值模型无法建立、结构表面振动信息难以获取。传统数值法(边界元法BEM、有限元法FEM、统计能量法SEA等)将受大网格计算量和复杂海洋环境和多物理场耦合等因素的严重限制，无法开展相关研究工作；理论解法一般采用镜像法原理，通过平整的镜像边界考虑浅海上下边界对结构声源的反射作用，但是这种方法的前提需要已知结构辐射声场的表达式，所以只能分析部分简单结构在简单波导下的声场问题，对复杂海洋波导下任意弹性结构的声辐射研究也无能为力。也有部分学者借鉴海洋声传播理论(如射线法、简正波理论、抛物线方程法和波数积分法以及混合方法等)的研究成果，类比分析结构声源在浅海波导下辐射声场的传播规律。虽然，随着距离的增加，具有几何尺寸的结构声源可近似等效为单个点源，其辐射声场的传播规律采用点源计算获取，这样的类比分析具有一定的借鉴价值。但是在近场直接忽略结构与流体、结构与环境的耦合作用，无法准确地进行弹性结构近场声场相关领域(声场测量、减振降噪和近场声辐射特性等)的研究。以上共同导致了无法从理论解和数值法的角度有效地开展波导下弹性结构声辐射研究，但它对海洋中水下结构声辐射、预报和识别有极为重要的意义，急需探索一种新的研究方法来解决海洋波导中结构声辐射问题。

发明内容

本发明提供了一种Pekeris波导下结构声振计算的FE/WSM方法，本发明对弹性结构和海洋波导环境适应能力强，计算准确高效，且易于在实际工程应用中推广，详见下文描述：

一种Pekeris波导下结构声振计算的FE/WSM方法，所述方法包括：

采用有限元法建立局域Pekeris波导下弹性结构声辐射多物理场耦合模型，进行三维Pekeris波导下弹性结构离散处理；通过有限差分法联立求解各个耦合模型，提取局部Pekeris波导下结构表面振动数据矩阵U；

采用波叠加法进行弹性结构内部虚拟点声源的空间优化布放，并基于球面波分解的精准镜像源法求解近场Green函数，获取各个虚拟源与结构表面离散节点的传递矩阵D，通过振动数据矩阵U和传递矩阵D的逆运算求解近场各个虚拟源的声源强度Q；

通过远场Green函数计算各个虚拟源到场点的单极矩阵T；将声源强度Q与对应的单极矩阵T相乘获取声场P；将各个离散虚拟源在远场场点产生的声场进行叠加求和，以等效为连续弹性结构的辐射声场。

进一步地，所述方法还包括：通过虚拟声源构成虚拟面，设置虚拟面与弹性结构的距离；

若弹性结构的表面为光滑刚性材料，距离为弹性结构的几何中心到表面距离的一半，且声源需在结构内部；

若弹性结构的表面为粘弹性材料，距离在最小结构波长之内。

其中，所述基于球面波分解的精准镜像源法求解近场Green函数具体为：

在Pekeris波导下边界液-液交互面上，通过波数域和空间域的坐标变换，将各个球面镜像源分解为各阶平面波和非均匀平面波；

将平面波和非均匀平面波与对应角度的界面反射系数相乘；采用自适应的分段积分牛顿-柯茨公式将各个角度对应下的反射声场积分，获取Green函数，进而求解近场声场。

其中，所述通过远场Green函数计算各个虚拟源到场点的单极矩阵T具体为：

根据频率与特征函数的变化关系，采用扰动理论进行复特征值的超越方程求解，获取远场Green函数，进而求解远场声场。

进一步地，所述方法采用有限元法建立局域Pekeris波导下弹性结构声辐射多物理场耦合模型具体为：

通过有限元法将声场方程、结构振动方程和特定边界条件联立求解，获取近场模型下结构振动和辐射声场的信息。

其中，所述方法基于球面波分解的精准镜像源法求解近场Green函数具体为：

采用自适应的分段积分牛顿-柯茨公式将各个角度对应下的反射声场进行积分，求解和计算近场Green函数。

本发明提供的技术方案的有益效果是：

(a)该发明对弹性结构形状和海洋波导环境类型适应能力强，采用有限元法建立海洋波导下弹性结构多边界耦合影响作用下的近场声辐射模型，因此快速准确计算任意复杂结构在不同海洋波导下的辐射声场，对结构几何类型、材料属性和复杂程度具有很强的适应性。然后，把局域下计算结构表面法向振速，作为波叠源强求解的输入条件，并结合近场Green函数求解各个虚拟源，最后再结合远场Green函数进行海洋波导下远场声场叠加计算。且在本发明可通过调整有限元建模、波叠加计算以及Green函数推导来适应不同的浅海环境，如本发明中，通过在Green函数求解过程中把声速和波数调整为复数来考虑液态海底对结构辐射声场的声吸收作用。进一步，若海底为质地较硬的弹性层或较为复杂的多孔材料层，则可通过只调整下连续边界，然后求解远近场Green函数，最后结合波叠加法求解在具有弹性或多孔海底的海洋波导下弹性结构辐射声场，且对于浅海楔形海底或深海波导下结构声辐射预报问题也同样使用类似处理方法，所以该方法可对不同海洋波导下的任意弹性结构辐射声场进行快速预报。

(b)该发明在进行计算时具有计算用时少、效率高的特点，因为由于海洋环境下辐射声场研究一般重点关注低频。所以有限元离散网格量少，可采用有限节点的法向振速信息表征结构表面振动特性，然后在内部布放有限个虚拟源，通过这些少量虚拟源产生声场进行叠加求和，得到结构等效的辐射声场。所以在相同计算硬件条件下，可以快速的计算传统数值法(有限元法、边界元法、统计能量法等)无法模拟或因计算量大而无法计算的复杂声辐射问题，且计算结果准确高效。

(c)该发明操作过程简单使用方便，实施工作量小，易于在理论研究和实际工程中推广。在实际应用中，其关键一步为只需在结构表面固定少量的加速度计，测量各个测点的法向振速，便可计算整个弹性结构的辐射声场，以实时监测水下弹性结构辐射声场。

与以往计算方法相比，该发明对弹性结构和海洋波导适应能力强，且计算结果准确效率高，实际操作简单易于推广。有效地解决目前在研究海洋波导下弹性结构声辐射预报时所遇到的计算量大、物理场众多和波导环境复杂等诸多瓶颈性问题。

附图说明

图1为本发明采用的Pekeris波导下弹性结构辐射声场FE/WSM预报模型；

图2(a)为多物理场耦合有限元数值法示意图；

图2(b)为弹性结构与周围流体的流固耦合；

图3为基于球面分解的镜像法示意图；

图4为本发明实施流程图；

图5为Pekeris波导下典型点源声传播示意图；

图6(a)为三种方法求解Green函数在频率为50Hz下声场对比；

图6(b)为三种方法求解Green函数在频率为100Hz下声场对比；

图6(c)为三种方法求解Green函数在频率为200Hz下声场对比；

图6(d)为三种方法求解Green函数在频率为400Hz下声场对比；

图7(a)为Pekeris波导下弹性圆柱壳声辐射示意图；

图7(b)为Pekeris波导下弹性圆柱壳声辐射有限元网格图；

图8(a)为频率为50Hz下圆柱壳表面振速分布云图；

图8(b)为频率为200Hz下圆柱壳表面振速分布云图；

图9(a)为频率为50Hz下，本发明声场计算结果与有限元计算结果对比图；

图9(b)为频率为200Hz下，本发明声场计算结果与有限元计算结果对比图。

表1为Pekeris波导环境参数配置表。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。

为了解决背景技术中存在的问题，有效地开展浅海波导下结构声辐射计算方法的研究，本发明实施例需要突破以往单一弹性结构振动和海洋声传播理论研究的技术壁垒，通过结构声学与海洋声学相结合，准确求解浅海波导下结构-流体-激励-边界综合系统中多个耦合子系统，快速进行海洋环境下结构辐射声场的传播计算。

实施例1

参见图1至图5，本发明实施例提供了一种基于FE/WSM的Pekeris波导下弹性结构声振计算方法，该方法包括以下步骤：

步骤一：采用有限元法建立局域Pekeris波导下弹性结构声辐射多物理场耦合模型，进行三维Pekeris波导下弹性结构离散处理；通过有限差分法联立求解各个耦合模型，提取局部Pekeris波导下结构表面振动数据矩阵U；

其中，在弹性结构表面与外部流体接触的耦合面上，通过结构表面法向的振动速度与外部流体介质的振动速度相同，建立起结构振动方程和流体域声学波动方程的相互作用关系(即流固耦合方程)。在浅海上边界，满足声压为零的边界耦合条件；液态海底下边界满足为声压和法向振速连续的耦合条件。无穷远处满足Smerfield辐射条件即远处辐射声压为零。

步骤二：采用波叠加法进行弹性结构内部虚拟点声源的空间优化布放，并采用基于球面波分解的精准镜像源法求解近场Green函数，获取各个虚拟源与结构表面离散节点的传递矩阵D，通过振动数据矩阵U和传递矩阵D的逆运算求解近场各个虚拟源的声源强度Q；

步骤三：在简正波理论求解远场Green函数过程中，采用扰动理论进行复特征值的求解，通过该远场Green函数计算各个虚拟源到场点的单极矩阵T；将各个虚拟源的声源强度Q与对应的单极矩阵T相乘，得到各个虚拟声源在场点的声场P；将各个离散虚拟源在远场场点产生的声场进行叠加求和，以等效为连续弹性结构的辐射声场。

本发明还可以包括：

1、在进行结构内部虚拟声源优化布放时，需要合理设置虚拟声源构成虚拟面与结构表面的距离：当结构表面为光滑刚性材料时，该距离为结构几何中心到结构表面距离的一半；当结构为粘弹性材料时，该距离在最小结构波长(结构纵波波速除以最大分析频率)之内。

2、不同于以往直接采用简单镜像源求解的Green，即直接采用球面波与平面波反射系数乘积再做叠加运算。本发明在采用镜像源法进行Green函数求解时，在Pekeris波导下边界液-液交互面上，通过波数域和空间域的坐标变换，把各个球面镜像源分解为各阶平面波和非均匀平面波；然后将这些平面波和非均匀平面波与对应角度的界面反射系数相乘；最后，采用自适应的分段积分牛顿-柯茨公式将各个角度对应下的反射声场积分。获取精准的Green函数，进而求解近场声场问题。

3、由于在该技术整个计算步骤中(即有限元数值模型、虚拟源强求解和声场计算)，均考虑了声场传播时介质的衰减吸收，导致采用简正波求解Green函数时出现了复平面求解多值超越方程。本发明实施例中根据频率与特征函数的变化关系，采用扰动理论进行复特征值的超越方程求解，得到较为准确且高效的Green函数，以进行远场声场快速计算。

实施例2

下面结合具体的实例对实施例1中的方案进行进一步地介绍，详见下文描述：

如图1所示，为本发明实施例的Pekeris波导下任意弹性结构低频声辐射的FE/WSM理论模型示意图。该理论模型由近场有限元法用于离散结构和远场波叠加用于计算声场构成。其计算步骤分主要包括：

1、Pekeris波导下波叠加法

如图1所示，流体密度为ρ_a，流体声速为c_a。弹性结构表面记为S，n为结构表面外法向矢量。P为外部流体中场点声压。Ω为结构内部与结构共形的等效源面，Q_i为Ω上声源强度，Q_i和P间的距离矢量为r_i，r_si和r_0i分别为结构中心到结构表面一点和虚拟源的距离。

在理论上，波叠加法与Helmholtz积分公式等效，场点的声压P可表示沿着虚拟源面Ω的积分：

其中，j是复数的虚数单位，ω为角频率，G(r_si-r_oi)是浅海波导中声场Green函数，通过Green函数考虑了浅海波导上下边界对结构和声场的作用。

在浅海环境下，Green函数需要满足以下方程：

其中，δ(r_si-r_oi)为狄利克雷函数，c(z)为介质声速。

求解以上Helmholtz方程便可得到对应浅海下的Green函数，不同的声传播模型对上述方程的数学处理不一样，另外在所用的方程的形式也不同。每一种声场计算模型都有其一定的适用条件，对于本发明实施例所涉及的Pekeris波导下结构声辐射计算问题，将采用镜像法和简正波法进行Green函数的求解，以高效准确求解远、近场声场问题。

通过线性Euler公式，将式(1)中的声压P表示为法向振速U具体为：

其中，Q(r_oi)为内部各个点源的声源强度。

经过离散化处理，将上式改写成矩阵形式，复杂结构源表面上某点的法向振速可由N个虚拟声源构建：

式中，

为求导算子，N为内部点源个数，Q_i(r_oi)为内部第i个点源的源强，u_n(r_s)为结构表面法向振速。

将上式改写成矩阵形式，表示为：

{U}_M×1＝{D}_M×ND{Q}_N×1 (5)

其中，M、N分别为结构表面上速度和虚拟源面点源的离散点数。

由于结构表面法向速度矩阵U是已知的，可通过FEM数值法获取；D为源点和场点的传递函数，一旦配置好等效源，传递函数

可通过相应浅海信道环境下的Green函数得到。通过已知U、D便可计算虚拟源的源强矩阵Q：[Q]＝[D]^-1[U](6)

其中，[D]^-1为[D]的Moore-Penrose广义逆。

获取虚拟源源强后，便可求解任意场点的声压：

P＝TQ＝T[D]^-1[U] (7)

其中，T为单极矩阵，定义为T＝jρωG(|r-r₀|)，_ρ为介质密度，r为结构几何中心到场点的距离，r₀为结构几何中心到点源的距离。

结构辐射声场一般关注于结构辐射声压级，定义为：

其中，p_ref＝1×10^-6Pa为水中的参考声压。

2、多物理场局域FEM数值计算

从信道下弹性结构辐射声场波叠加计算式(6-7)可知，该方法关键的第一步是获取结构表面振速，且其获取的准确性将决定整个声场计算过程的精度。不同于其他流体环境下结构声振问题，浅海信道下结构表面的激励条件和声场将受上下界面反射声的重要影响，形成了激励力-结构-流体声场-浅海边界相互耦合作用的复杂系统，且由于海底边界类型多、声学参数复杂、地形多变等因素限制，加大了进行浅海信道下结构表面声振信息获取的难度和工作量，解析法难以推导，而试验法成本高。由于浅海下结构表面振动特性只与距离范围内的边界反射声相关，且针对信道下低频声辐射问题由于被动吸声技术发展和高频衰减快，有限元网格量相对较少，所以采用有限元法建立多边界下结构近场低频声辐射局域数值模型获取结构表面振动信息是可行的。

1)流固耦合方程

建立如图2(a)所示的浅海信道下弹性结构声辐射局域数值计算模型，其中Ω_a为浅海流体介质域，四周Ω_p为声学无反射层即完美匹配层(Perfectly Matched Layer，简称PML)域，其宽度为h_PML；PML内外边界为Γ_I、Γ_O，Γ为流固耦合边界，Γ_l、Γ_h分别为声场与浅海海面与海底交互作用边界。根据对应边界下的连续条件，建立结构与流体，结构与边界耦合方程。

对Helmholtz方程进行权重积分并结合高斯理论，可写出声学有限元方程为：

(K_a+jωC_a-ω²M_a){p_i}＝{F_i}on S (9)

其中，{F_i}为声学激励，M、K和C分别为质量矩阵、刚度矩阵和阻尼矩阵，下标a为声学系统。

类似声学有限元方程推导，对于弹性结构，其有限元振动方程为：

(K_s+jωC_s-ω²M_s){d_i}＝{F_si}in Ω_a (10)

其中，M_s、K_s和C_s分别为结构网格上没有受到约束(位移d_i)部分的刚度矩阵、质量矩阵和阻尼矩阵；{F_si}为结构上的激励载荷。

在结构表面与外部流体接触的耦合面上，如图2(b)所示，满足的边界条件为结构表面法向的振动速度与外部流体介质的振动速度相同，可写出结构与流体的耦合方程为：

其中，刚度矩阵K和阻尼矩阵C、质量矩阵M均为n×n阶矩阵。下标a、s和c为声学系统、力学结构和耦合系统；定义耦合矩阵K_c、M_c为

且

T为矩阵转置运算，N网格数量。n_se为结构与流体接触的结构网格数量，{n^e}为结构网格的法向向量；ω＝2πf为角频率，f为频率(Hz)，ρ₀为海水密度；u_i、p_i为位移和声压幅值，F_st、F_at分别为结构、流体介质声的耦合激励载荷。

2)Pekeris波导上下边界耦合连续条件

海洋波导的海面边界通常为Dirichlet边界，满足的边界条件为界面声压为零：

p_a(x,y,z)|_z＝0＝0 (12)

其中，为了区分海水和海底，下标a和b分别标示海水流体层和海底层。

对于液态海底，满足的边界条件为声压p(x,y,z)连续，法向振速v(x,y,z)连续

p_a(x,y,z)＝p_b(x,y,z) (13)

v_an(x,y,z)＝v_bn(x,y,z) (14)

其中，v_an和v_bn分别为耦合面上、下层介质的速度。

海洋波导四周边界为无限远边界，有限元数值法采用PML技术进行模拟，PML通过对控制方程增加吸收系数转换为吸收层的控制方程，为了简化方程描述，令x轴为x₁轴、y轴为x₂轴，利用分离变量可写出频域下的PML方程：

其中，σ_i为吸收系数，v_i，p_i为匹配层域的速度和声压幅值。

3)四周无限远边界

采用PML对边界进行处理后，使在边界上满足Smerfield远场熄灭条件：

p(x,y,z)|_r＝∞＝0 (16)

在边界上通过声吸收使声压为零，达到边界没有反射声以模拟波导四周的无限大空间。Pekeris波导下受激弹性结构的振动与声辐射机理为：作用于结构表面的激励力能量传遍整个弹性结构，导致结构因为应变产生弯曲振动产生激励的简谐波，这种振动通过式(11)所示的流固耦合方程导致周围流体介质产生压缩和伸张运动，引起介质声场的传播。且声场传播需要满足特定的浅海边界条件，例如在上边界需要满足式(12)所示的Dirichlet边界，下边界满足式(13-14)所示的液态连续边界，然后四周边界为无限远边界，采用PML模拟，满足式(16)所示的Smerfield远场熄灭条件。通过有限元法把声场方程、结构振动方程和特定边界条件联立求解，计算获取近场模型下结构振动和辐射声场的信息。

3、Pekeris波导下Green函数

从大陆架到大陆架坡的过程中，海底岩基表面覆盖一层非凝固态物质，这些物质由大量江河流入的沉积物，所以把海底沉积物模拟为液态是合理的，且低频辐射声场入射海底会与海底沉积物产生不可忽略的能量损耗，这里必须引入合理的衰减效应，以提高点源声传播(即Green函数)求解精度。

浅海信道下点源声传输函数即Green函数的求解方法有多种(2式)，如射线理论、简正波方法、波数积分方法和抛物方程等方法。根据各种方法的适用条件，且在波叠加源强求解涉及近场声场问题(4式)，远场声场叠加计算涉及远场问题(5式)，为了提高波叠加整个计算精度，这里采用虚源法和简正波法分别进行近、远场的声场计算。

镜像法是按射线理论而来的，如图3所示，理想浅海下点源产生的声场是各个声线进行干涉叠加的总和，即：

其中，V₁,

分别为海面海底的平面波发射系数，理想海面边界通常为软边界即，V₁＝-1，Pekeris波导下海底并非理想的硬海底，而是液态海底，海底反射系数并非保持不变，与式与角度相关的，

为球面波分解的各阶平面波与海底的入射角(并非球面波与海底法线的角度)；k_a为海水中声场波数，

z_m1＝2Hm+z_s-z，z_m2＝2H(m+1)-z_s-z，z_m3＝2Hm+z_s+z，z_m4＝2H(m+1)-z_s+z，H为水深。m为本体源与前三个虚源镜像次数，一般当m达到几十阶以后，场点总声场计算结果会很快稳定地于某一个值。

但是，不同于硬海底条件时的虚拟源镜像叠加，即界面反射系数始终为1。Pekeris波导下，液态海底的反射系数不仅与入射角度相关，还与波的形式相关。且对平面波的声场反射系数更为熟悉，而式(17)中声场为各个虚拟源的球面波叠加。这里为了提高近场声场计算的准确性以提高波叠加法的精度，很有必要把各个球面波分解为各阶平面波和非均匀平面波。通过双重傅里叶积分，可把球面波分解为平面波，即：

其中，

为直角坐标系下平面波矢k_a与z轴的夹角，

为第一类汉克耳函数，θ_i为镜像源O_mn的传播路径R_mn与z轴的夹角。

上述积分的上下限为复平面无穷积分，采用常规方法难以准确求解，且求解效率低。本发明采用自适应的分段积分牛顿-柯茨公式把各个角度

对应下的反射声场进行积分，可准确求解和计算Green函数。

通过式(18)把球面波分解平面波，对于海底边界，反射系数为：

其中，

β＝c_a/c_b，

为平面波的入射角，ρ_a、ρ_b分别为海水与海底的密度，c_a、c_b分别为海水与海底的声速。

当海底考虑损耗效应时，海底边界反射系数不在为一实数，而变成以虚数，因为在当海底存在吸收时，海底声速为：

其中，α_b为沉积物衰减系数，与采用每波长衰减的吸收系数的关系为α_λ(dB/λ)＝27.3α_b。

采用简正波求解Pekeris波导下Green函数，把式(2)改写为深度方向和距离方向的柱坐标下，最终方程可以写成：

其中，z_s为源位置，z为场点位置，δ(r)和δ(z-z_s)分别为距离和深度方向上的狄拉克函数。

这里，把Pekeris波导设为均匀水体，即海水密度和声速为常数，则上式改为：

其中，k_a为水中声场波数。

采用分离变量求解以上非齐次方程的解，将G(r,z)＝Φ(r)Z(z)带入方程(22)，然后采用完备正交基处理后，便可将声场描述为各个简正模式的叠加，最后声压可写为：

其中，r、z分别为轴对称坐标下水平方向和竖直方向，z_s为声源位置，z为场点位置，j为复数虚部，

为汉克函数，ξ_n为水平方向特征值；Z_n(k_zz)解形式由海面、海底边界决定，k_z为竖直方向的特征值。

Z_n(k_zz)定义为：

Z_n(k_zz)＝B sin(k_zz)+Ccos(k_zz) (24)

其中，B、C分别为正弦和余弦函数的幅度值。

海面通常设为Dirichlet边界条件，满足的边界条件为：

Z_n(k_zz)|_z＝0＝0 (25)

在半无限液态海底模型交互面上，满足声学边界方程为：

其中，g(ξ_n)由液态海底声学特性决定，通过g(ζ_n)也可建立与界面反射系数的关系，即为：

同样，当海底考虑声吸收时，可以通过复声速得出海底的复波数即：

将式(24)分别带入式(25-26)，便可得出特征方程：

求解式(29)的特征值，为复数域的超越方程。常规的数值方法难以准确求解。这里，根据频率与特征函数的变化关系，本发明采用扰动理论进行复特征值的超越方程求解，得到较为准确且高效的Green函数，以进行远场声场快速计算。

由理论模型可知，如图(4)所示。Pekeris波导下结构声辐射波叠加法计算过程主要分为三步：第一步，采用多物理场耦合有限元理论建立近场局部声辐射模型，联立求解流固耦合方程[式(11)]、上下边界和无穷远耦合条件[式(12-15)]，计算获取结构表面的声振信息；第二步，采用式(6)和近场Green函数[式(17)]进行源强Q的求解；第三步，利用式(7)和远场Green函数[式(23)]计算Pekeris波导下结构的辐射声场P。

实施例3

算例1：Pekeris浅海波导下点源声传播计算

为了验证通过镜像法和简正波法求解Green函数的准确性，图5所示，为建立的Pekeris波导下点源声传播模型示意图，模型中结构和Pekeris波导环境参数如表1所示。点源深度为z_s＝10m,场点深度为z＝15m，相邻间距为1m。波导上边界为压力释放边界，下边界为含有声学吸收的半无限液态空间(密度，声速，吸收系数)，波导深度为30m，海水密度声速分别为。分别采用镜像法(image source method，ISM)求解的近场Green函数即式(17)、简正波法(normal mode，NM)推导的远场Green函数即式(23)计算了该环境下点源的辐射声场声压级，并与相同条件下的波数积分法(wavenumber integration，WI)计算的精准解进行了对比分析。

表1

如图6(a)-(d)所示，分别为采用镜像法和简正波法计算了50、100、200和400Hz的声场声压级，并与相同条件下精准的波数积分解(采用成熟的Scooter程序计算)进行了对比，计算距离范围为0-2000m。可看出，由于球面波在平整面上的反射会产生平面波和非均匀平面波(泄漏模式)，而在简正波解中未考虑这些高阶泄漏模式，是一种近似解，其近场声场计算精度与精准WI解有一定的差异，不能由于近场声场问题的准确求解。但泄漏模式会按照距离出现1/r²衰减，且在2-3波长之内基本完全衰减。所以达到一定距离以后，简正波解也能取得很高的计算精度，且计算高效，所以很适合处理远场声场计算问题。而镜像法通过把上下边界设为镜像边界，通过把多个镜像反射声源与本征声源产生的声场进行叠加，且对于非理想边界，通过把球面波分解为各阶平面波，然后通过平面波反射系数处理了液态海底边界对声场的反射作用。所以该镜像法是一种类似波数积分法的精确解，在整个计算距离上的声场计算精度与波数积分法解完全吻合，且求解过程较为简单相对于波数积分法，可以很好地处理近场源强求解。在源强求解过程中，由于点源与结构表面距离在几个波长以内，属于近场声学问题，而在浅海结构辐射声场波叠加声场计算中，通常关注于远场辐射声场传播的规律，为了同时兼顾FE/WSM的计算精度和效率，采用镜像法和简正波法进行波叠加过程中源强求解和声场计算。

算例2：Pekeris浅海波导下弹性结构声辐射计算

为了显示该方法对于三维结构声辐射计算的优势，进行了三维浅海波导下非轴对称的圆柱壳声辐射验证。如图7(a)所示，为Pekeris波导下大型弹性圆柱壳声辐射有限元模型，浅海环境、参数与上述相同，圆柱壳长为10m、半径1m、厚度0.01m，中心深度为15m，材料也为4340型钢(密度ρ_s＝7850kg/m³,杨氏模量E_s＝2×10¹¹Pa,泊松比u_s＝0.33)，Pekeris波导水深为H＝30m，环境参数、分析频率和场点选取与上述验证模型1相同，在圆柱壳上表面轴线中心施加垂直向下的简谐点激励力，激励力幅值为1N。在有限元计算中，如图7(b)所示，流体域网格划分大小满足网格大小小于六分之一波长，弹性结构以及结构附近的流体域的局部模型的网格划分更为精细，在波导四周和海底底部采用六层网格的PML层来吸收声波，已达到界面处无反射声，模拟四周无限远边界和海底无限大半空间域，通过合理的有限元网格划分，提高有限元计算的准确度和计算效率。同时，在FE/WSM计算过程中，表面振速也是通过局部模型计算获取，近场局部模型的有限元计算精度要求更为严格，网格划分质量要求更高。

通过有限元近场模型计算获取弹性圆柱壳表面S上的振速分布云图，如图8(a)和8(b)所示，分别为50Hz和200Hz对应的圆柱壳表面振速分布，然后均匀选取了表面176个节点的法向振速，结构内部虚拟源面Γ均匀分布160个点源，且与结构表面共形，几何尺寸比为2:1(即半径0.5m，轴线长度5m)，因为圆柱壳结构为三维非轴对称条件下的辐射声源，所以需要更多的虚拟源以达到理想的计算精度。目前，尚无高效可靠的Pekeris波导下弹性圆柱壳声辐射计算的准确方法，本发明采用FE/WSM计算了三维Pekeris波导下弹性圆柱壳辐射声场，并与相同条件下有限元计算结果进行了对比，由于为三维声场问题，为了减小网格计算量，计算距离设为0-500m，频率为50Hz和200Hz。

如图9(a)和9(b)所示，分布为50Hz和200Hz下，本发明计算结果与有限元计算计算结果对比图。可看出，除了在近距离处与有限元解存在一些偏差，因为非对称弹性结构振动和辐射声场较为复杂，通过有限点源来刻画结构的声振特性确实会现在一定的精度问题，但整体计算精度是满足计算要求的.

而且，该计算算例为贴近于实际尺寸的大型弹性结构声辐射模型，以往常规方法无法准确计算该模型。但是本发明方法能够准确快速的计算各个方位和频率下的辐射声场。该算例中，对200Hz辐射声场计算用时为24.91分钟(计算硬件参数为戴尔工作站，型号Precision 5820Tower，主频为the Intel(R)W-2133 3.60GHz,内存为96GB)，局部近场有限元离散时产生的计算自由度为2.11×10⁶。相反，如果采用有限元计算该三维声辐射问题时，其200Hz的网格计算量为1.33×10⁷，计算用时约为13小时。如果，当加大计算距离或升高计算频率，将在常规有限元计算中导致计算自由度急剧增加，计算机内存溢出而无法计算，但是该发明通过局部有限元建模，网格计算量少，且通过有限个虚拟源结合高效的Green函数，便可快速计算任意距离和频率下的辐射声场。所以，该发明不仅能够准确计算三维弹性结构声辐射，且计算量独立于计算距离和频率，计算代价低、效率高。

本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制，只要能完成上述功能的器件均可。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种Pekeris波导下结构声振计算的FE/WSM方法，其特征在于，所述方法包括：

采用有限元法建立局域Pekeris波导下弹性结构声辐射多物理场耦合模型，进行三维Pekeris波导下弹性结构离散处理；通过有限差分法联立求解各个耦合模型，提取局部Pekeris波导下结构表面振动数据矩阵U；

采用波叠加法进行弹性结构内部虚拟点声源的空间优化布放，并基于球面波分解的精准镜像源法求解近场Green函数，获取各个虚拟源与结构表面离散节点的传递矩阵D，通过振动数据矩阵U和传递矩阵D的逆运算求解近场各个虚拟源的声源强度Q；

通过远场Green函数计算各个虚拟源到场点的单极矩阵T；将声源强度Q与对应的单极矩阵T相乘获取声场P；将各个离散虚拟源在远场场点产生的声场进行叠加求和，以等效为连续弹性结构的辐射声场。

2.根据权利要求1所述的一种Pekeris波导下结构声振计算的FE/WSM方法，其特征在于，所述方法还包括：通过虚拟声源构成虚拟面，设置虚拟面与弹性结构的距离；

若弹性结构的表面为光滑刚性材料，距离为弹性结构的几何中心到表面距离的一半；

若弹性结构的表面为粘弹性材料，距离在最小结构波长之内，且所有虚拟声源在结构内。

3.根据权利要求1所述的一种Pekeris波导下结构声振计算的FE/WSM方法，其特征在于，所述基于球面波分解的精准镜像源法求解近场Green函数具体为：

在Pekeris波导下边界液-液交互面上，通过波数域和空间域的坐标变换，将各个球面镜像源分解为各阶平面波和非均匀平面波；

将平面波和非均匀平面波与对应角度的界面反射系数相乘；采用自适应的分段积分牛顿-柯茨公式将各个角度对应下的反射声场积分，获取Green函数，进而求解近场声场。

4.根据权利要求1所述的一种Pekeris波导下结构声振计算的FE/WSM方法，其特征在于，所述通过远场Green函数计算各个虚拟源到场点的单极矩阵T具体为：

根据频率与特征函数的变化关系，采用扰动理论进行复特征值的超越方程求解，获取远场Green函数，进而求解远场声场。

5.根据权利要求1所述的一种Pekeris波导下结构声振计算的FE/WSM方法，其特征在于，所述方法采用有限元法建立局域Pekeris波导下弹性结构声辐射多物理场耦合模型具体为：

通过有限元法将声场方程、结构振动方程和特定边界条件联立求解，获取近场模型下结构振动和辐射声场的信息。

6.根据权利要求1所述的一种Pekeris波导下结构声振计算的FE/WSM方法，其特征在于，所述方法基于球面波分解的精准镜像源法求解近场Green函数具体为：

采用自适应的分段积分牛顿-柯茨公式将各个角度对应下的反射声场进行积分，求解和计算近场Green函数。

Images