一种预报水下声辐射的无网格弱式径向点插值方法

吴绍维，柯 磊，韩国文

(1. 重庆交通大学 航运与船舶工程学院，重庆 400074；2. 船舶动力工程技术交通运输行业重点实验室，武汉 430063)

水下声辐射问题是近年来非常活跃的研究领域，基于网格的方法是现阶段计算水下声辐射的主要手段，如边界元法(boundary element method, BEM)和有限元法(finite element method, FEM)。BEM和FEM均可直接求解内部声学问题，对于外部声问题，BEM在无穷远处满足Sommerfeld辐射条件，是求解无限域声学问题的有效手段[1]。但BEM在处理大型问题时，由于相关系数矩阵为非对称和非稀疏矩阵，导致计算效率变低[2]。FEM是求解大型问题非常有效的方法[3]，然而标准的FEM不能直接用于求解无限域中的外部声学问题，需对问题域进行处理，常采用人工边界截断无限域获取有限计算域，为确保声波向外传播过程中的自由衰减[4]，需在人工边界处施加无反射边界条件，如DtN(Dirichlet-to-Neumann)边界条件、完美匹配层、无限元和吸收边界条件等。在标准的FEM中，存在数值色散误差[5-6]，细化网格不能有效降低该误差[7]。近年来提出了一些基于光滑数值和无网格的方法来软化有限元模型“过硬”的刚度[8-10]。借助这些方法，色散误差得以减小。

用无穷傅里叶级数表示的DtN边界条件是一种精确的非局部边界条件[11-12]，解析构建了Dirichlet变量与Neumann变量之间的关系。在计算时，需将无穷级数截断为有限项数N，当N不足够大时，无法确保解的唯一性和可解性[13-14]。因此，截断形式的DtN边界条件不再精确，虽然通过增加项数可解决该问题，但对高频或大尺度人工边界计算量过大。为克服该缺陷，MDtN(modified Dirichlet-to-Neumann)边界条件被提出，对于任意的N值，可确保声场求解的唯一性和可解性[15]。与截断的DtN边界条件相比，MDtN边界条件使用少量的项数可获得更高的精度。此外，在声场计算中无法预先确定具体的项数以达到期望的精度。因此，MDtN边界条件更具优势。

与基于网格的方法相比，新兴的无网格法无需使用节点连接信息或网格进行变量插值或近似[16]，已被广泛应用于处理一些具有挑战性的数值和工程问题。无网格点插值法是一种全局弱式法[17]，该方法采用多项式点插值法(polynomial point interpolation method, PPIM)构造形函数，使之具有Kronecker-delta函数性质，可以像有FEM一样直接施加本质边界条件。然而，PPIM形成的矩矩阵可能是奇异的。使用两阶矩阵三角化算法可消除该奇异问题[18]，但由于PPIM形函数的不相容性，使得该方法对不规则分布节点的鲁棒性较差，且局部支持域中含有过多节点会形成过高阶的插值多项式，导致PPIM形函数剧烈振荡。为消除奇异性，基于径向基函数(radial basis function，RBF)的径向点插值法(radial point interpolation method, RPIM)被提出[19]，其对任意分布的节点具有良好的稳定性和鲁棒性，具有比FEM更快的收敛速度。

RPIM已被用于内部声学问题，能获得比FEM更精确的结果[20]。作为一种基于域离散的方法，RPIM需与其他数值方法联合处理无限域声问题[21]。在水下声辐射计算领域，为抑制中频数值色散误差，基于光滑有限元耦合DtN边界的方法被提出，包括混合光滑有限元法、稳定点基光滑有限元法[22]及边基梯度光滑法[23]。这些耦合方法有效提高了水下中频段声辐射和声散射计算精度。在光滑有限元框架下，Xu等[24]提出了单元基光滑径向基点插值耦合DtN边界条件方法计算水下声辐射，将RPIM引入水下噪声计算，与传统的FEM相比，该方法在计算精度、收敛速度和抑制数值色散误差方面具有优势。基于RPIM和MDtN边界条件的各自优点，本文提出了一种基于RPIM和MDtN边界条件相耦合的无网格弱式方法来精确预报水下声辐射。研究了无网格插值中的尺寸参数和形参数对声场预报精度的影响规律，确定了可用取值范围。结果表明，与FEM相比，所提方法在计算精度、收敛速度以及计算效率方面具有优势，对声波数的敏感度显著降低。

1 水下声辐射计算理论背景

1.1 控制方程

在无限大自由场V中，稳态简谐结构振动辐射声波的声压p(x)满足如下Helmholtz方程

∇2p(x)+k2p(x)=0,x∈V

(1)

式中： ∇2为拉普拉斯算子；k=ω/c为声波波数，ω为结构振动圆频率，c为流体介质中的声传播速度。结构边界可分为Dirichlet型边界ΓD和Neumann型边界ΓN，ΓD和ΓN满足ΓD∪ΓN=Γ及ΓD∩ΓN=Ø，其定义为

p=pD, onΓD

(2)

∇p(x)·nb=-ikρcv(x), onΓN

(3)

式中：p=pD为ΓD上的已知声压；v(x)为ΓN上给定的法向振速；nb为边界外法向；ρ为流体介质密度； i2=-1。

对exp(iωt)形式的简谐时间因子，为确保声波向外传播过程中的自由衰减和声场求解的唯一性要求，需满足如下Sommerfeld辐射条件[25]

(4)

式中：r为场点到坐标原点的距离；β为维度参数，对二维和三维问题分别取1和2。若p(x,t)=p(x)exp(-iωt)，则∂p/∂r-ikp=o(r-β)。通常只用第二项表示Sommerfeld 辐射条件，在无穷远处，此辐射条件与ρc阻抗条件等效，即∇p·n∞+ikp=0。

1.2 MDtN边界条件

为获取有限计算域，采用人工边界ΓR(对2D和3D问题分别为半径为R的圆和球)截断无限问题域，有限计算域以Γ为内边界，ΓR为外边界，如图1所示。在人工边界上施加DtN边界条件∂p(x)/∂n=-Mp(x)，式中M为DtN映射算子，其作用是在ΓR上建立Dirichlet量声压p与Neumann量∂p/∂n的关系。从而，最初的声场求解问题可等效为如下两个声场求解问题。

(5)

(6)

图1 使用一组离散分布场点表示的有限计算域Fig.1 Finite computational domain represented by a set of discrete distributed field nodes

为完整起见，简要介绍DtN边界条件理论，更为详细的理论和推导可参考Keller等的研究。对二维声问题，采用极坐标(r,θ)可解析求解出式(6)对应问题的声压为

(7)

(8)

其中，

(9)

mn(θ-θ′)=

(10)

在实际计算中，需将式(8)中的无穷级数截断为N个有限项，则式(8)变为

(11)

式中，MN为截断后的DtN映射算子。从而截断后的DtN边界条件不再精确，当N取值较小时，无法确保问题的可解性和数值解的唯一性。为克服该缺陷并提高计算精度，MDtN边界条件被提出，以确保N取任意值时具有唯一的解。通过在式(8)右边加上和减去Bp(R,θ)，同时将Bp(R,θ)求和并截断，得到如下MDtN边界条件

(12)

(13)

则MDtN边界条件可确定为

(14)

其中，

(15)

2 RPIM-MDtN计算理论

2.1 形函数构造

采用ph(x)表示声压p(x)的试函数，W(x)表示检验函数，式(1)两边乘以W(x)，并在VI上积分得到

(16)

根据格林第一公式，得到如下方程

(17)

式中：n为包裹VI的边界外法向；在Γ上，n=-nb。将式(17)代入式(16)得到

(18)

将式(3)和式(12)代入式(18)得到如下弱式方程

(19)

用散布在VI及其边界上的一组离散场点表示该有限计算域及其边界，VI中任意计算点x处的声压可使用场点声压值插值确定。

(20)

式中：Ri(x)为RBF；pj(x)为基于空间坐标x=[x,y,z]T表示的单项式；n和m分别为RBF基函数和多项式的项数；a=[a1,a2,…,an]T和b=[b1,b2, …,bm]T为未知的常数系数向量；RT=[R1(x),R2(x),…,Rn(x)]，pT=[p1(x),p2(x), …,pm(x)]。

为求解式(20)中的未知系数a和b，通过在计算点x的支持域Ωs中的n个场点处满足式(20)来构造n个方程，并结合如下m个附加方程

(21)

得到如下系统方程

(22)

式中：ps=[p1,p2, …,pn]T为声压值向量；

(23)

(24)

对任意分布的场点，R-1通常存在，且式(20)通常使用低阶多项式，则RPIM中不会出现奇异问题。求解式(22)得到系数a和b，将其代入式(20)得到

(25)

式中，NT(x)=[N1(x),N2(x), …,Nn(x)]为n个场点的RPIM形函数。在RPIM形函数的构造中，需确定局部支持域，通常采用影响域来确定计算点x的支持域中包含的场点，若计算点位于某场点的影响域中，则该场点参与构造该计算点的形函数。在后续部分，为了简化起见，采用圆形影响域。根据建议，采用rw=αsd确定影响域半径rw，式中αs表示场点xi的影响域的无量纲尺寸，该参数对计算精度至关重要，对于给定类型的问题，通常采用数值试验来确定αs取值。

2.2 RPIM-MDtN离散系统方程

在加权残差公式中，采用RPIM构造的形函数N(x)作为检验函数，则式(19)可表示为

(26)

经离散化处理，建立如下矩阵形式的离散系统方程

[K-k2M+Kb]{p}=F

(27)

(28)

(29)

(30)

(31)

式中： {p}为声压值的向量；Cn(x)=[cos(nθ),sin(nθ)]。以上这些矩阵是基于场点i和j的全局索引组装，计算这些系统矩阵需用到一组全局背景单元，其不依赖于场点，不参与变量插值，对插值精度无影响。只有当场点i和j同时属于至少一个积分点的支持域时，Kij≠0和Mij≠0，否则，Kij=0且Mij=0。因此，全局系数矩阵K和M为稀疏矩阵。

3 数值算例

采用数值算例对所提方法进行验证。通过使用具有解析解的浸没无限长圆柱数值试验确定用于水下声辐射计算的影响域无量纲尺寸和形状参数的取值，并对该方法的精度、收敛性以及计算效率开展研究。进一步通过浸没无限长圆柱局部声辐射模型、舵形结构和二维潜艇结构的声辐射算例，对所提方法进行检验。

3.1 浸没周向简谐振动无限长圆柱声辐射

无限长圆柱半径为a=0.25 m，声介质为水，密度取ρ=1 000 kg/m3，声速为c=1 500 m/s，以坐标原点为中心构建R=0.5 m圆形人工边界，如图2所示。在圆柱的湿表面上(r=a)指定周向Dirichlet边界条件p(θ)=cos(nθ)，对应的归一化解析解为

(32)

采用εreal和εimag分别表示声压实部和虚部数值均方根误差

(33)

(34)

图2 浸没周向简谐振动无限长圆柱声辐射模型Fig.2 Submerged infinitely long cylinder with circumferential oscillation

3.1.1 影响域无量纲尺寸对精度的影响

对一类基准问题，通常需进行数值试验来预先确定αs的可用取值，本节通过分析声压预报误差随αs的变化来研究αs对精度的影响规律。使用一组平均节点间距为d=0.023 m的不规则分布场点离散有限计算域，采用496个二次四边形单元作为全局背景单元实施积分运算(如图3所示)，构造形函数时在式(20)添加线性多项式。声波波数满足ka=π/2，使得a等于一个波长。为了确保无网格全局弱式法的积分精度，对于二维问题，积分点nq与场点nf的数量应满足nq>2nf/3。由于nf=1 171，为精确计算系统矩阵系数，在背景单元的每个方向上使用5个Gauss-Legendre积分点。图4和图5分别为第二阶(n=1)和第三阶(n=2)辐射模式下人工边界处的εreal和εimag随αs的变化曲线，其中ac=0.4和q=1.03(MQ)，ac=2(EXP)，η=3(TPS)，N=7(MDtN)。

图3 使用任意分布的场点表示的计算域及用于积分的背景单元Fig.3 Representing computational domain using arbitrarily distributed nodes and performing integration with the use of background cells

图4 第二阶辐射模式下εreal和εimag随αs的变化Fig.4 εreal and εimag against αs for second circumferential mode

图5 第三阶辐射模式下εreal和εimag随αs的变化Fig.5 εreal and εimag against αs for third circumferential mode

由图中曲线可见，声压预报误差随αs改变而变化，αs取值对所提方法的精度具有显著的影响。为了获得精确的计算结果，需合理选取αs。当影响域尺寸过小时(αs≤1.5)，矩阵G为奇异矩阵，导致计算失败。因此，图中没有给出对应αs≤1.5的计算结果。采用EXP-RBF作为基函数时，用所提方法获得的结果在2≤αs≤6范围可接受；当αs<2或αs>6时，精度降低。对于MQ-RBF和TPS-RBF，在αs≥2范围内，所提方法可以产生精确的结果，但随着αs的增大，并不能有效提高计算精度，其原因可能是αs与形状参数不匹配，对于该问题还需作进一步研究。然而，较大的影响域将显著增加计算时间，尤其是对于大规模问题。因此，需保持计算精度和计算耗时之间的平衡。结果还表明，相比EXP-RBF，采用MQ-TPF和TPS-RBF作为基函数具有更高的精度。为了不增加计算耗时并获得较高的计算精度，在后续算例中，对EXP-RBF取αs=3，对MQ-RBF和TPS-RBF取αs=4.5。

3.1.2 形状参数对精度的影响

为研究形状参数对所提方法精度的影响，图6～图9给出了人工边界处的εreal和εimag随形参数的变化。由图中曲线可知，计算精度受形状参数的影响，需合理选取形状参数的值以获得精确的计算结果。当采用MQ-RBF时，由图6可知，αc可取值范围为0≤αc≤5；图7中的曲线表明，形状参数q对计算精度具有显著影响，q=j(j=1,2,3,…)时，RPIM中的G矩阵为奇异矩阵，导致错误的计算结果。对于EXP-RBF，由图8可见，αc在0.07～2.50内可产生较精确的计算结果。当使用TPS-RBF时，从图9给出的结果可见，η应在1≤η≤9范围取值，需注意，当η取值接近2j(j=1,2,3,…)时，矩阵G的条件数变大，导致计算精度急剧变差。根据以上结果，在后续研究中，基函数的形状参数取值为：对MQ-RBF，取αc=0.5和q=1.1；对EXP-RBF，取αc=1.9；对TPS-RBF，取η=3.5。

图6 对应MQ-RBF的εreal和εimag随αc的变化Fig.6 εreal and εimag against αc corresponding to MQ-RBF

图7 对应MQ-RBF的εreal和εimag随q的变化Fig.7 εreal and εimag against q corresponding to MQ-RBF

图8 对应EXP-RBF的εreal和εimag随αc的变化Fig.8 εreal and εimag against αc corresponding to EXP-RBF

图9 对应TPS-RBF的εreal和εimag随η的变化Fig.9 εreal and εimag against q corresponding to TPS-RBF

3.1.3 有限项数N的影响

为研究MDtN边界条件中的截断项数对所提方法精度的影响规律，图10给出了在第五阶辐射模式下人工边界处的误差随有限项数N的变化，同时还给出了使用DtN边界条件产生的数值误差，其中波数k满足ka=2π。结果表明：与截断的DtN边界条件相比，使用MDtN边界条件可获得更高的计算精度，尤其是对于较小的N取值，计算精度提升明显。对于截断的DtN边界条件，通过设定N≥kR，虽能够克服适定性问题，但由于通常无法预先确定所需项数来达到期望的精度，计算结果仍可能严重偏离真实值。与截断的DtN边界条件不同，MDtN边界条件对积分算子内核中包含的低阶模态是精确的，而局部算子包含了与高阶模态对应的近似非反射边界条件。因此，MDtN边界条件对应的误差在N=5时快速收敛到稳定值，在后续算例中，取N=6。

图10 MDtN和DtN边界条件对应的误差随N的变化Fig.10 Errors obtained from MDtN and DtN maps versus N

3.1.4 收敛性和计算效率研究

为研究收敛性，采用四种平均场点间距(d=0.011, 0.020, 0.032,0.043)进行数值试验。图11给出在第三阶辐射模式下所提方法和FEM的收敛曲线，其中k满足ka=2π，两种方法采用相同的场点配置，使用二次有限元单元和和对应的节点。图11中的曲线表明，在相同的配置下，与FEM相比，所提出的方法具有更高的收敛率，能够产生更精确的结果；当使用MQ-EBF或EXP-RBF时，收敛过程不稳定，一个可能的原因是无量纲尺寸与形状参数不匹配。到目前为止，还没有可用的方法从理论上确定无量纲尺寸和形状参数的取值，只能通过数值试验获得，需进一步从理论上确定其最优值。

在相同配置下，通过对比所提方法与有限元法的计算耗时来开展计算效率研究，所有方案及程序均在相同的计算环境中运行(Intel®Xeon®W-2255 CPU 3.7 GHz和RAM 64 GB)。在效率研究中，采用TPS-RBF作为基函数。图12(a)给出了在ka=2π时，CPU耗时随d的变化。结果表明：在相同配置下，所提方法需更多的计算时间。其原因是，对所提方法，必须为所有积分点支持域内的场点构造形函数，而在FEM中，已预先确定覆盖积分点的单元的FEM形函数；此外，与FEM相比，所提方法需使用更多的场点信息来组装系统矩阵，导致系统矩阵的带宽较FEM大。众所周知，有效的数值方法应以短的计算耗时获取高计算精度，因此将计算耗时与精度结合能合理评计算效率，图12(b)给出了CPU计算时间随声压实部误差的变化。由结果可见，当εreal=17.5%时，该方法的计算耗时约为FEM的39%，从而当要求高精度时，所提方法具有计算效率优势。

图11 所提方法与FEM的收敛性比较Fig.11 Comparison of convergence curves between the proposed scheme and FEM scheme

图12 所提方法与FEM的计算效率比较Fig.12 Comparison of computational efficiency between the proposed scheme and FEM scheme

3.2 浸没无限长圆柱局部声辐射

为进一步检验采用所提方法预报声场的准确性，采用浸没无限长圆柱的局部声辐射模式进行数值试验，对于该算例，在圆柱边界局部范围(-φ≤θ≤φ,φ=π/6)指定p=1，在其余边界处p=0，如图13所示。几何尺寸、相关参数和平均场点间距与上述算例相同，对应的归一化解析解为

(35)

图14和图15分别为k=4π和k=8π，声压沿人工边界的分布。由结果可见，数值解与解析解总体上吻合良好，高声波数对应的声压值与解析解略有偏离，但计算结果仍可接受，与EXP-RBF相比，使用MQ-RBF和TPS-RBF作为基函数可获得更精确的结果。

图13 浸没无限长圆柱局部声辐射Fig.13 Non-uniform acoustic radiation by a sector of cylinder

图14 当k=4π时人工边界处的声压分布Fig.14 Sound pressure on artificial boundary when k=4π

图15 当k=8π时人工边界处的声压分布Fig.15 Sound pressure on artificial boundary when k=8π

3.3 舵形结构声辐射

浸没于水中的刚性舵形结构的形状和尺寸如图16所示，水介质密度ρ=1 000 kg/m3，声速c=1 500 m/s，在结构的湿表面施加vn=10-5m/s的Neumann边界条件，以半径为R=4 m的圆作为人工边界获取有限计算域，并采用1 193个双线性四节点单元对其进行离散(含1 277个节点)。为进行对比，在相同场点配置下，给出FEM计算结果。由于该问题不存在解析解，采用细化网格(68 677个四节点单元，69 357个节点)对应的FEM结果作为参考值。图17和图18和分别给出kL=4π和kL=8π时声压沿人工边界的分布，其中RPIM采用TPS-RBF作为基函数进行插值。由图中结果可见，在给定频率处，RPIM-MDtN获得的声压分布与参考值吻合良好。对于较小的波数，FEM-MDtN计算结果与参考值具有良好的一致性，但随着声波数的增加，在kL=8π时，FEM预报的声压严重偏离参考值。结果表明，所提方法对声波数敏感度低，在水下声辐射计算中具有更高的精度。

图16 刚性舵形结构的几何形状与尺寸Fig.16 Shape and geometry of rigid rudder-shaped structure

图17 当kL=4π时人工边界处的声压分布Fig.17 Sound pressure on artificial boundary when kL=4π

图18 当kL=8π时人工边界处的声压分布Fig.18 Sound pressure on artificial boundary when kL=8π

3.4 潜艇结构声辐射

本节采用浸没于水中的二维潜艇结构的声辐射试验来检验所提方法的实用性，水的密度ρ=1 000 kg/m3，声速c=1 500 m/s，潜艇形状和尺寸如图19所示。对于潜艇，需评估其辐射噪声水平，敌人声纳通常距离被探测物体较远，在所提RPIM-MDtN方法中，在确定人工边界处的声压值后利用式(7)可计算任意远场点处的声压。在实际预报潜艇噪声时，潜艇的艉部和艏部船体的振动模态取决于分析频率[26]，为简单起见，仅考虑位于潜艇尾部的动力装置引起的船尾振动。在此算例中，在尾部施加如图19所示的振动速度为vn=10-5m/s的Neumann边界条件。以原点为中心，半径为R=50 m的圆作为人工边界截断无限问题域，采用10 545个双线性四节点单元(10 841个节点)离散计算域，并在人工边界处选取观测点(θ=π/2)。由于该声辐射问题不存在解析解，使用细化的网格(177 820个四节点单元，179 005个节点)对应的FEM结果作为参考值。

图19 二维海底形结构的几何形状和尺寸Fig.19 Size and shape of submarine-shaped structure

图20～图21给出了f=100 Hz及400 Hz时，声压沿人工边界的分布，其中采用TPS-RBF构造形函数。为进行对比，在相同配置下给出FEM的计算结果(10 545个双线性四节点单元，10 841个节点)。由图中结果可见，在低频f=100 Hz时，RPIM-MDtN和FEM-MDtN方案均可产生良好的结果，随着频率的增加，RPIM-MDtN预报的声压依然非常接近参考值，而FEM-MDtN预报的声压在f=400 Hz时明显偏离了参考值。为进一步检验所提方法的精度，研究观测点处的10～400 Hz的声压频率响应，其中频率间隔为1.0 Hz。RPIM-MDtN和FEM-MDtN预报的声压及参考值绘于图22。由结果可见，相比FEM-MDtN，RPIM-MDtN在全频率范围内具有更高的精度且在高频处预报的结果与参考值仍然具有良好的一致性，而FEM-MDtN在高频处预报的结果严重偏离参考值。结果再次表明：与FEM相比，所提方法在水下声辐射预报中能够获得更高的计算精度，对声波数的敏感度显著降低，能够有效抑制中频色散误差。

图20 当f=100 Hz时人工边界上的声压分布Fig.20 Sound pressure on artificial boundary at f=100 Hz

图21 当f=400 Hz时人工边界上的声压分布Fig.21 Sound pressure on artificial boundary at f=400 Hz

图22 观测点处的声压频率响应Fig.22 Sound pressure frequency response at observation point

4 结 论

基于径向点插值法和修正的Dirichlet-to-Neumann边界条件提出了一种预报水下结构声辐射的无网格弱式方法，通过算例对该方法进行了研究，主要结论总结如下：

(1) RPIM-MDtN能够精确且稳定地计算无限域水下声辐射，无需使用网格或节点连接信息来进行场变量插值。

(2) 与传统的FEM相比，RPIM-MDtN能够获得更精确的声场预报结果，在收敛性方面表现得更好，对声波数的敏感度显著降低，能够有效抑制中频色散误差。

(3) 由于PRIM形状函数构造，在相同场场点配置下，RPIM-MDtN方法比标准的FEM需要更多的计算时间，当要求高计算精度时，所提方法具有计算效率优势。

(4) 相比EXP-RBF，MQ-RBF和TPS-RBF具有更高的精度；TPS-RBF具有稳定的收敛性，而EXP-RBF和MQ-RBF的收敛过程不稳定，对于该问题还需作进一步的研究。

(5) 在所提方法中，影响域无量纲尺寸和形参数取值均采用数值试验获取，可能存在无量纲尺寸与形状参数不匹配的问题，需从理论上进一步研究其优化取值。