三维浅海下圆柱壳声辐射预报方法研究

三维浅海下圆柱壳声辐射预报方法研究钱治文; 何元安; 商德江; 肖妍; 史博【期刊名称】《《船舶力学》》【年(卷),期】2019(023)010【总页数】11页(P1266-1276)【关键词】浅海波导; 弹性圆柱壳; 耦合计算; 波叠加法; 传输函数; 声辐射特性【作者】钱治文; 何元安; 商德江; 肖妍; 史博【作者单位】天津大学海洋科学与技术学院天津300072; 中国船舶工业系统工程研究院北京100036; 哈尔滨工程大学水声工程学院哈尔滨150001【正文语种】中文【中图分类】TB5320 引言我国海洋领土大多属于水深小于200 m 的典型浅海环境，潜器在浅海所处的声场环境并非理想的自由场，存在水面和水底的反射、散射作用[1]，其声场特性与结构本身和浅海环境两者密切相关。

研究浅海下结构辐射声场特性，对开展浅海环境下结构振动辐射噪声实时监测预报、声学隐身性能评价以及浅海环境参数对辐射声场影响分析等研究领域具有重要的理论和应用价值，是今后我国水声技术领域长期关注的热点和难点问题之一。

然而，目前关于水下弹性结构耦合振动与声辐射问题的研究，大多考虑无界或半空间流体域下的近距离辐射声场问题[2-3]。

对于浅海环境下结构辐射声场问题，国内外很少涉及该方面的研究，因该研究在试验法、解析法和数值法研究过程中均遇到了难题，有效开展工作较为困难；试验法周期长成本高，数值法受网格划分以及计算量的严重制约，解析解法目前只能针对理想浅海环境下简单结构声场问题[4-5]。

为了解决浅海下弹性结构声辐射研究所遇到的这些难题，本文引入了一种对结构和流体环境适应性很强且计算高效的研究方法:波叠加法（Wave Superposition Method，简称WSM）。

Koop-mann 等[6]最先提出了利用波叠加法进行结构辐射声场的计算，Miller 和Fahnline 等随之对稳定性、计算精度和效率进行了分析[7-8]。

波叠加法以适应性强、计算效率高的优势逐步地成为了研究弹性结构声辐射的重要方法之一，近年来受到了国内学者极大的关注，并广泛应用于自由场或半空间流体环境下的近场声全息[9-11]、结构声辐射预报[12-13]和噪声源识别[14-15]等研究领域。

但采用波叠加法进行浅海下结构声辐射的研究却很少，因为浅海下结构辐射声场问题涉及流体、结构和浅海上下边界等多个物理场之间的相互耦合作用，导致结构表面振动速度难以准确获取，从而无法进行波叠加法的准确计算。

也有部分研究把浅海下结构辐射源视为点声源，但这样直接忽略了结构与流体、结构与环境的耦合作用和近场声辐射特性。

导致目前尚无高效可靠的研究方法进行浅海环境下弹性结构声辐射预报的有效研究，但它对浅海下结构声辐射预报、源识别和噪声控制等具有重要的意义，急需探索一种新方法来解决浅海下结构声辐射预报问题[16]。

本文提出了一种三维浅海波导下弹性结构辐射声场预报的理论方法。

该方法首先利用有限元法建立了浅海波导下结构声辐射多物理场耦合有限元模型；计算获取在考虑流体环境和浅海波导上下边界耦合影响下结构表面振动数据，然后采用波叠加法和浅海声场传输函数进行虚拟源源强求解，最后进行弹性结构辐射声场的等效叠加计算。

经有限元法和解析解法验证该方法的正确性后，研究了典型弹性圆柱壳在理想浅海波导下的辐射声场特性，为今后开展复杂浅海环境下大型弹性结构声辐射预报研究提供了一种新途径。

1 理论模型1.1 浅海波导下点源声传播模型考虑声源为水下简谐点源，设海水为理想介质，压力p 满足非齐次亥姆霍兹方程式中:▽2为拉普拉斯算符，c 为流体声速，r→为声源到场点的矢量距离，δ r→（）为狄拉克函数，A 为点源强度幅值，j 为复数虚部，ω 为角频率。

由于声场具有圆柱对称性，选用柱坐标下，以通过点源且垂直向下的方向为z 轴，垂直于z 轴的方向为r 轴，方程可以写成[1]式中:k0 为流体介质中声波波数，定义为k0=ω/c，z0为声源在z 轴上的位置。

采用简正波理论求解式，可得浅海波导下点源声场表达式p （r, z ）即浅海波导传输函数G （r, z ）为式中:ρ1为海水密度，z 为场点位置；为汉克函数，ξn 为水平方向特征值，kz 为竖直方向特征值，满足为格林函数，且为一个正交函数族，其解由海面、海底声学边界共同决定。

海面通常为Dirichlet 边界条件，满足的边界条件为海底是一个反射和散射边界，典型海底边界模型有Neumann 边界、Rayleigh 定律、Sommerfeld 模型、Cauchy 边界和地声模型等。

在Neumann 边界、Rayleigh 定律和Sommerfeld 模型的交互面上，满足声学边界方程为式中:g（ξ ）由各浅海海底的声学特性决定，定义为，m1=ρ1/ρ2，ρ2为海底密度，ξ 为海底水平方向波数，k2为海底波数，V 为海底界面反射系数。

对于复杂的Cauchy 边界和地声模型，其反射系数求解可参照文献[17]，然后便可求解声场传输函数。

1.2 波叠加法模型建立如图1 所示的三维浅海下结构声辐射波叠加法计算理论模型，流体密度为ρ1，声速为c1；结构中心和点源在海面的投影为o 和o′，o-xy 和o′-x′y′与海面重合，o-xyz 为全局坐标系，o′-x′y′z′为各虚拟源的局部坐标系，oo′连线延长线为y′轴，o′y′与oy 的夹角为θ；振动体表面记为S′，n 为结构表面的外法向矢量，zs 为结构中心在竖直方向的位置；S 为与结构面共形的等效源面，Q 为虚拟源，z0i 为各虚拟源在局部坐标中的位置；P 为场点位置，Q 和P 间的距离记为rp。

图1 三维浅海波导下波叠加法原理示意图Fig.1 Principle of wave superposition method in three-dimensional shallow water waveguide波叠加法与Helmholtz 积分公式等效，经过离散化处理后，全局坐标系下结构表面某一点的法向速度un 可由N 个简单声源构建[18]式中:qi是虚拟面上第i 区域上源强强度，rs为结构表面一点空间位置矢量，r0i 为虚拟面各点的空间位置矢量，G）为三维全局坐标下的声场传输函数。

把（7）式写成矩阵形式，可表示为结构表面法向速度矩阵U 已知，可通过有限元计算或实验获取，D 为虚拟源和场点的传递矩阵，定义为一旦配置好等效源，便可计算虚拟源的源强矩阵Q。

其中，［D］-1 为［D ］的广义逆矩阵。

为了求解传递矩阵的广义逆［D］-1，需要对其进行奇异值分解式中:S、V 分别为M 阶、N 阶方阵，Σ=diag （σ1, σ2, σ3…σn ）为包含n 个非奇异值σi 的对角阵。

通过奇异值分解求逆矩阵获取辐射体结构内部虚拟源源强后，便可计算任意一点的声压式中:T 为单极矩阵，为三维全局坐标下声场传输函数。

从（9）式和（11）式可知，传递矩阵D 和单极矩阵T 中涉及的传递函数G （x, y,z ）均为全局坐标下三维声场传递函数，而1.1 节求解的声场传递函数为局部柱对称坐标下二维声场传递函数G （r, z ），需要进行如下空间坐标转化。

浅海下各虚拟源产生的空间声场满足轴对称性，三维局部坐标系o′-x′y′z′下声场可降为二维轴对称局部坐标系o′-rz′下声场，根据局部二维坐标下各点源到场点水平距离r 和竖直距离z′求解各点源对场点的声场作用G （r, z′ ）即G （x′ , y′, z′ ）（其中x′=rcosϑ，y′=rsinϑ，ϑ为r 轴与x 轴夹角），然后通过局部坐标系与全局坐标系的转换关系建立各点源在全局坐标系下的声场作用G （x, y,z ）。

首先将局部坐标系声场函数G （x′ , y′, z′ ）沿x 方向、y 方向和z 方向分别平移Δx、Δy 和Δz，得到全局坐标系下声场G （xi, yi, zi ），平移变换矩阵为然后再将全局坐标系下声场函数G （xi, yi, zi）绕z 轴顺时针旋转角度θ，得到与全局坐标系o-xyz下的声场函数G （x, y,z ），旋转变换矩阵为通过局部坐标o′-rz′到全局坐标o-xyz 的转换，可求得各局部坐标下点源声场函数G （r, z′ ）对全局坐标系下结构表面的作用G （x, y,z ）。

由（9）式和（11）式可知，该方法计算的基本思路为:首先通过有限元法建立浅海波导下弹性结构声辐射数值模型，计算获取结构表面法向速度U，然后通过波叠加法以及传输函数G 进行虚拟源源强Q 求解，最后结合单极矩阵T 便可计算在任意场点的辐射声场P。

2 数值分析2.1 方法准确性与高效性分析如图2 所示，建立了浅海波导下脉动球有限元模型，波导深度h=50 m，海面为Dirichlet 边界，海底为Neumann 边界，流体四周无限大边界。

小球中心深度为zs=25 m，半径为r0=3 m，表面施加均匀振速ua，数值计算获取表面振动数据，然后采用波叠加法计算场点声场信息。

根据脉动球与点源的等效关系为则脉动球壳在浅海波导下的声场等效为浅海波导下点源声场，浅海波导下点源在任意场点的声压表达式为[1]式中:为简正波阶次n 的最大正整数为简正波特征值。

图2 浅海波导下脉动球声辐射有限元模型Fig.2 The finite element model of spherical shell in shallow water waveguide如图3 所示，利用本文浅海波导下辐射声场波叠加计算式即（11）式计算了浅海波导下脉动球辐射声压级随距离的变化曲线（各场点深度为h2=30 m），并分别与解析解法即（14）式、有限元法计算结果进行了对比分析，验证了该方法计算声场的可靠性。

图3 波叠加法计算结果验证Fig.3 The verification of wave superposition method图4 FFP 与Kraken 计算结果对比Fig.4 The comparison between FFP andKraken虽然采用简正波法存在一定的近场声场计算精度问题，但从远场波叠加法声场计算结果验证来看，其计算精度是满足要求的。

为了说明采用简正波法的声场波叠加计算精度问题，本节采用最为精准的波数积分法计算程序FFP 与简正波计算程序Kraken 进行了计算结果对比分析。

计算环境以Pekeris波导为例，海水密度ρw=1 024 kg/m3，声速cw=1 500 m/s；半无限液态海底密度ρw=1 800kg/m3，海底声速cw=2 500 m/s；海水深度h=200 m，声源深度zs=100 m，接收点深度zr=50 m，如图4 所示，分别采用FFP 和Kraken 计算了点源在浅海波导内的传播损失曲线对比图。