船舶力学

频率（Hz）	500	1000	1500	2000	2250
FEM-T4	2.95%	10.36%	16.40%	31.98%	49.15%
MRPIM	2.58%	0.65%	3.14%	8.90%	11.32%

频率（Hz）	500	1000	1500	2000	2250
FEM-T4	2.95%	10.36%	16.40%	31.98%	49.15%
MRPIM	2.58%	0.65%	3.14%	8.90%	11.32%

每波长节点数	12	10	8	6
FEM-T4	7.27%	16.53%	33.62%	49.15%
MRPIM	1.47%	3.11%	6.02%	11.32%

每波长节点数	12	10	8	6
FEM-T4	7.27%	16.53%	33.62%	49.15%
MRPIM	1.47%	3.11%	6.02%	11.32%

水下弹性目标声学分析的改进径向点插值无网格计算模型

PDF下载

游翔宇 ¹^,²^,³ , 姚宇 ¹^,³ , 李威 ⁴

船舶力学 | 流体与结构声学 2026,30(4): 648-658

收起

船舶力学 |流体与结构声学 2026 , 30 (4) : 648 -658

水下弹性目标声学分析的改进径向点插值无网格计算模型

全屏

游翔宇¹^,²^,³, 姚宇¹^,³, 李威⁴

作者信息

^1.长沙理工大学水利与海洋工程学院，长沙 410114

^2.上海交通大学海洋工程全国重点实验室，上海 200240

^3.水沙科学与水灾害防治湖南省重点实验室，长沙 410114

^4.华中科技大学船舶与海洋工程学院，武汉 430074

通讯作者:

姚　宇（1982–），男，博士，教授，通讯作者，E-mail：yaoyu821101@163.com

作者简介:

游翔宇（1992–），男，博士，讲师

李　威（1975–），男，博士，教授

A modified radial point interpolation meshless model for acoustic analysis of underwater elastic objects

Xiang-yu YOU¹^,²^,³, Yu YAO¹^,³, Wei LI⁴

Affiliations

^1.School of Hydraulic and Ocean Engineering, Changsha University of Science & Technology, Changsha 410114, China

^2.State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

^3.Key Laboratory of Water-Sediment Sciences and Water Disaster Prevention of Hunan Province, Changsha 410114, China

^4.School of Naval Architecture and Ocean Engineering, Huazhong University of Science and Technology, Wuhan 430074, China

出版时间: 2026-04-15 doi: 10.3969/j.issn.1007-7294.2026.04.013

文章导航

摘要

收起

经典有限元法（FEM）在较高计算频率下会受频散误差效应的影响，出现显著的精度劣化。此外，有限元解的可信度高度依赖于网格质量，尤其在具有陡峭梯度场或复杂几何特征的区域，高质量网格是获得可靠结果的必要前提，但此类网格的生成过程往往具有计算密集性和耗时性。针对上述问题，本研究提出一种改进径向基点插值无网格计算方法（MRPIM）用于三维声固耦合分析。该方法通过引入新型插值节点选择方案，在单个积分单元内构造了连续声压近似函数，从而有效降低了数值积分误差。结果表明，改进方法在计算精度上较原有方法（RPIM）具有显著提升，且有更优的计算效率；数值分析进一步揭示，该方法与有限元法的计算成本具有可比性，显示出良好的工程应用潜力。

关键词

有限元法 / 三维声学分析 / 无网格法 / 频散误差 / 新型插值节点选择方案

Abstract

收起

The classical finite element method (FEM) is known to suffer from significant accuracy degradation under large wavenumber conditions due to dispersion error effects. Furthermore, the reliability of FEM solutions heavily depends on mesh quality, particularly in regions with steep gradient fields or complex geometric features. While high-quality meshes are essential for obtaining credible results, their generation remains computationally intensive and time-consuming. To address these limitations, this study proposes a modified radial point interpolation meshless method (MRPIM) for three-dimensional acoustic-structure interaction analysis. By introducing a novel interpolation node selection scheme, the method constructs continuous approximation functions within individual integration cells, effectively reducing numerical integration errors. Results demonstrate that the modified method achieves notable improvements in computational accuracy compared to the conventional method (RPIM) while exhibiting superior efficiency. Additionally, the computational overhead of MRPIM is comparable to that of FEM, highlighting its promising potential for engineering applications.

Key words

finite element method / three-dimensional acoustic analysis / meshless method / dispersion error / novel interpolation node selection scheme

引用本文

游翔宇, 姚宇, 李威. 水下弹性目标声学分析的改进径向点插值无网格计算模型. 船舶力学, 2026 , 30 (4) : 648 -658 . DOI: 10.3969/j.issn.1007-7294.2026.04.013

Xiang-yu YOU, Yu YAO, Wei LI. A modified radial point interpolation meshless model for acoustic analysis of underwater elastic objects[J]. Journal of Ship Mechanics, 2026 , 30 (4) : 648 -658 . DOI: 10.3969/j.issn.1007-7294.2026.04.013

正文

收起

0　引　言

收起

水下目标探测技术中，声学探测是可靠的技术手段之一。通过分析目标的辐射或散射声波，可解读目标的有无、位置或材料属性等信息^[1−2]。为了分析简单形状弹性目标的声学特性，学者提出了许多物理声学方法，如分波序列法、Sommerfeld-Watson变换法和共振散射理论等。然而工程中普遍存在的是复杂形状目标，上述方法对此无能为力，因此研究数值离散方法变得十分必要。经典数值离散方法如有限差分法FDM（Finite Difference Method）、有限元法FEM（Finite Element Method）和边界元法BEM（Boundary Element Method）等在声学计算中被广泛使用，可处理各类复杂形状目标^[1]。但这些方法的可靠性高度取决于网格的数量和质量，导致工程分析中耗费大量人力及时间修整耦合界面或梯度剧烈变化区域的网格^[3]。因此，发展无网格的数值离散方法成为解决复杂工程应用的有效途径之一。

基于伽辽金弱形式的无网格方法最早诞生于上世纪90年代，如无单元伽辽金法（Element-Free Galerkin Method）^[4]和再生核质点法RKPM（Reproducing Kernel Particle Method）^[5]等。之后，Liu和Gu基于径向插值基函数，发展了一种单变量的径向点插值无网格方法RPIM（Radial Point Interpolation Method）^[6]。该方法具备上述两种方法缺失的Kronecker Delta函数性质，因而能直接处理本质（声压）边界条件；并且，该方法的形函数实际上只有与空间距离有关的一个变量，在分析三维问题时更加简单灵活。鉴于上述优点，国内外有学者将RPIM应用于求解内场^[4,7−8]和外场声学问题^[9−10]，同等离散尺度下相比，FEM取得了更精确的波动数值解。然而RPIM生成的声压近似函数不连续，理论上可能导致数值不稳定现象，因此有学者结合广义梯度光滑技术发展了分区光滑径向点插值无网格法以解决此问题^[8]。但广义梯度光滑技术使用了常值光滑函数（即光滑域面积），会抵消形函数的高阶性质，导致精度难以进一步提高。

为提高RPIM的精度及稳定性，作者提出基于积分网格中心点的插值节点选择方案^[11]，构造改进径向点插值无网格法MRPIM（Modified RPIM），用于求解二维声学问题，结果表明，该方案可有效提高原有RPIM法的计算精度及效率。为进一步解决实际声学工程问题，在本文中作者基于改进方法构造三维无网格声固耦合计算模型，并对其计算精度、效率和收敛性等进行了研究，结果表明，相比于经典FEM和RPIM，MRPIM各方面的表现更优，具备良好的实用性。

1　改进径向点插值无网格法

收起

径向点插值无网格法采用离散场节点进行插值，但仍需借助背景积分网格对控制方程的弱形式进行全局积分。如图1所示，与FEM不同，无网格法离散后的场节点（蓝色圆圈）可不与三角形积分网格（黑色虚线）的顶点重合，意味着无网格法可动态调整节点位置而无须重构网格，离散时更加灵活。此外，传统插值节点选择方案中基于积分网格的高斯积分点（x形）构造形函数，表明一个积分网格需要多次进行插值节点选取和“矩”矩阵（由径向基函数和多项式基函数在支撑域内所有节点上的值构成的一个系数矩阵）求逆操作，计算成本较高，且形成的近似函数不连续，可能导致数值积分过程的不稳定；而改进的插值节点选取方案基于积分网格的中心点（三角形）构造形函数，一个积分网格仅需一次选取插值节点，计算成本大大降低。

径向点插值无网格法使用了多二次（Multiquadrics）径向基函数进行插值，其表达式如下

$ {R}_{i}\left(r\right)={\left[{r}^{2}+{\left({\alpha }_{{\mathrm{c}}}{d}_{{\mathrm{c}}}\right)}^{2}\right]}^{q} $

(1)

式中：

$ r=\sqrt{{\left(x-{x}_i\right)}^{2}+{\left(y-{y}_i\right)}^{2}+{\left(z-{z}_i\right)}^{2}} $

表示空间两点之间的欧氏距离，d_c为节点平均间距，α_c和q为影响计算精度的待定参数。基于后续的数值试验，确定α_c = 0和q = 0.95可在三维声学分析中得到较高的计算精度。

如图1所示，在积分点的支撑域（彩色实线方形，半边长r_s = 3d_c）中，选取域内的所有节点进行插值

$ p^h(\boldsymbol{x})=\boldsymbol{R}(\boldsymbol{x}) \boldsymbol{a}+\boldsymbol{P}(\boldsymbol{x}) \boldsymbol{b}=\sum_{i=1}^n R_{i}(\boldsymbol{x}) a_i+\sum_{j=1}^3 P_j(\boldsymbol{x}) b_j $

(2)

式中：P_j = [1 x y z]^T表示为通过分片测试而补充的线性多项式基函数，a和b为对应基函数的广义坐标，n为支撑域内所有节点的数量。

将所有节点的坐标及节点声压代入式（2）中，可得

$ \boldsymbol{p}={\boldsymbol{R}}_{\text{M}}\boldsymbol{a}+{\boldsymbol{P}}_{\text{M}}\boldsymbol{b} $

(3)

式中：

$ \boldsymbol{R}_{\mathrm{M}}=\left[\begin{array}{cccc}R_1\left(r_1\right) & R_2\left(r_1\right) & \cdots & R_n\left(r_1\right) \\R_1\left(r_2\right) & R_2\left(r_2\right) & \cdots & R_n\left(r_2\right) \\\vdots & \vdots & \ddots & \vdots \\R_1\left(r_n\right) & R_2\left(r_n\right) & \cdots & R_n\left(r_n\right)\end{array}\right] $

(4)

$ \boldsymbol{P}_{\mathrm{M}}=\left[\begin{array}{cccc}1 & x_1 & y_1 & z_1 \\1 & x_2 & y_2 & z_2 \\\vdots & \vdots & \vdots & \vdots \\1 & x_n & y_n & z_n\end{array}\right] $

(5)

为保证插值结果的唯一性，额外的约束条件如下

$ \boldsymbol{P}_{\text{M}}^{\text{T}}\boldsymbol{b}=\boldsymbol{0} $

(6)

联立式（3）和式（6），并对“矩”矩阵求逆，可得

$ p^h(\boldsymbol{x})=\left[\begin{array}{ll}\boldsymbol{R}(\boldsymbol{x}) & \boldsymbol{P}(\boldsymbol{x})\end{array}\right] \underbrace{\left[\begin{array}{cc}\boldsymbol{R}_{\mathrm{M}} & \boldsymbol{P}_{\mathrm{M}} \\\boldsymbol{P}_{\mathrm{M}}^{\mathrm{T}} & \boldsymbol{0}\end{array}\right]^{-1}}_{\text {“矩”矩阵的逆}}\left[\begin{array}{l}\boldsymbol{p} \\\boldsymbol{0}\end{array}\right]=\left[\begin{array}{ll}\boldsymbol{\varPsi}(\boldsymbol{x}) & \boldsymbol{M}(\boldsymbol{x})\end{array}\right]\left[\begin{array}{l}\boldsymbol{p} \\\boldsymbol{0}\end{array}\right] $

(7)

式中：

$ \boldsymbol{\varPsi }\left(\boldsymbol{x}\right)=\boldsymbol{R}\left(\boldsymbol{x}\right)\boldsymbol{R}_{\text{M}}^{-1}+\left(\boldsymbol{P}\left(\boldsymbol{x}\right)-\boldsymbol{R}\left(\boldsymbol{x}\right)\boldsymbol{R}_{\text{M}}^{-1}{\boldsymbol{P}}_{\text{M}}\right){\boldsymbol{S}}^{-1}\boldsymbol{P}_{\text{M}}^{\text{T}}\boldsymbol{R}_{\text{M}}^{-1} $

为径向基形函数，

$ \boldsymbol{M}\left(\boldsymbol{x}\right)=\left(\boldsymbol{R}\left(\boldsymbol{x}\right)\boldsymbol{R}_{\text{M}}^{-1}{\boldsymbol{P}}_{\text{M}}-\boldsymbol{P}\left(\boldsymbol{x}\right)\right){\boldsymbol{S}}^{-1} $

为多项式基形函数（不参与节点插值），

$ \boldsymbol{S}=\boldsymbol{P}_{\text{M}}^{\text{T}}\boldsymbol{R}_{\text{M}}^{-1}{\boldsymbol{P}}_{\text{M}} $

。由上述插值过程可知，径向基形函数保留了Kronecker Delta函数性质，因此能直接处理声压边界条件。

在MRPIM中，由于使用了改进的插值节点选择方案，式（7）中“矩”矩阵的逆在每个积分网格仅需计算一次。而使用传统的插值节点选择方案时，“矩”矩阵的逆在每个积分网格的高斯积分点处计算一次。MRPIM作为高阶数值方法，需采用高阶高斯积分法才能保证积分精度。在本文中使用了16点四阶精度的高斯积分法在四面体网格中对弱形式进行积分，使用改进插值节点选择方案之后，可将求逆次数降低至传统方案的1/16，将显著提高形函数的计算效率。不仅如此，传统方案中每个积分点处选取不同的插值节点进行插值后，声压近似函数在积分网格内是不连续的，会导致额外的数值积分误差。作为对比，改进方案中在不同积分点处使用了相同的一套插值节点进行插值，保证了积分网格内声压近似函数的连续性，从而可相对提高计算精度。

2　水下弹性目标声学控制方程及完美匹配层技术

收起

对于线性稳态声学问题，亥姆霍兹方程是其控制方程，如下式

$ \Delta p+{k}^{2}p=0 $

(8)

式中：Δ是拉普拉斯算子，k = ω/c表示声波波数，ω为圆频率，c为媒介的声波波速，p为声压波动。

将式（8）乘加权函数w，并对整个问题域Ω进行积分，再运用格林散度定理可得弱形式

$ \int_{\varOmega} (\nabla w \cdot \nabla p-k^2 w p )\mathrm{~d} \varOmega-\int_{\varGamma} w(\nabla p \cdot \boldsymbol{n}) \mathrm{d} \varGamma=0 $

(9)

式中：

$ \varGamma $

为问题域的Neumann边界，

$ \boldsymbol{n} $

为边界

$ \varGamma $

的外法向单位矢量。

为了吸收外传声波，满足Sommerfeld远场辐射条件，本文引入了基于复坐标映射的完美匹配层技术PML（Perfectly Matched Layer）。如图2所示，将无限计算域截断为有限计算域，并在有限域外围设置PML层，α为层厚度。

各坐标轴方向上的复坐标映射关系如下

$ \tilde{x}_i=T\left(x_i\right)=\left\{\begin{aligned}&x_i-{i} \int_\eta^0 \sigma(\eta) \mathrm{d} \eta, \\&\qquad \text { 当 }-R_i-\alpha \leqslant x_i< -R_i \\&x_i, \quad \text { 当 }-R_i \leqslant x_i \leqslant R_i \\&x_i-{i} \int_0^\eta \sigma(\eta) \mathrm{d} \eta, \\&\qquad \text { 当 } R_i< x_i \leqslant R_i+\alpha\end{aligned}\right. $

(10)

式中：i表示坐标轴下标x、y和z，相对坐标

$ \eta =\left| {x}_i\right| -{R}_i $

，阻尼函数

$ \sigma \left(\eta \right)=7{\left(\eta /\alpha \right)}^{\text{2}} $

，具体配置及参数选取依据见文献[12]。

引入PML后，式（9）相应变换为

$ \int_{\varOmega_{\mathrm{p}}} \nabla w \cdot \nabla p-k^2 w p \mathrm{~d} \varOmega+\int_{\tilde{\varOmega}_{\mathrm{PML}}} \tilde{\nabla} w \cdot \tilde{\nabla} p-k^2 w p \mathrm{~d} \varOmega-\int_{\varGamma} w(\nabla p \cdot \boldsymbol{n}) \mathrm{d} \varGamma=0 $

(11)

式中：

$ \tilde{\nabla } $

表示相对于复坐标的梯度算子。

代入RPIM形函数

$ \boldsymbol{\varPsi }_\text{f}^{} $

后，式（11）可简写为

$ \left(\boldsymbol{K}_{\mathrm{f}}-k^2 \boldsymbol{M}_{\mathrm{f}}\right) \boldsymbol{p}=\boldsymbol{F}_{\mathrm{f}} $

(12)

式中：

$ {\boldsymbol{K}}_{{\mathrm{f}}}={\left(\nabla \boldsymbol{\varPsi }_{{\mathrm{f}}}^{}\right)}^{\text{T}}\nabla \boldsymbol{\varPsi }_{{\mathrm{f}}}^{} $

表示声流体刚度矩阵，

$ {\boldsymbol{M}}_{{\mathrm{f}}}=\boldsymbol{\varPsi }_{{\mathrm{f}}}^{\text{T}}\boldsymbol{\varPsi }_{{\mathrm{f}}}^{} $

表示声流体质量矩阵，

$\boldsymbol{F}_{\mathrm{f}}=\displaystyle\int_{\varGamma} \boldsymbol{\varPsi}_{\mathrm{f}}^{\mathrm{T}}\left(\nabla \boldsymbol{\varPsi}_{\mathrm{f}} \boldsymbol{p} \cdot \boldsymbol{n}\right) \mathrm{d} \varGamma $

表示外部声激励。

当考虑弹性体的声学特性时，还需引入弹性体的简谐振动响应方程^[1]

$ \left({\boldsymbol{K}}_{{\mathrm{e}}}-{\omega }^{2}{\boldsymbol{M}}_{{\mathrm{e}}}\right){\boldsymbol{u}}_{{\mathrm{e}}}={\boldsymbol{F}}_{{\mathrm{e}}} $

(13)

式中：K_e表示弹性体刚度矩阵，M_e表示弹性体质量矩阵，F_e表示外部力激励。

问题域内总声压p可表示为入射声压p_i与散射（或辐射）声压p_s之和，如下

$ p={p}_{{\mathrm{i}}}+{p}_{{\mathrm{s}}} $

(14)

在声固耦合面上，弹性体表面位移u_e应与声流体质点位移u_f相等

$ {\boldsymbol{u}}_{{\mathrm{e}}}{\boldsymbol{n}}_{{\mathrm{f}}}={\boldsymbol{u}}_\text{f}{\boldsymbol{n}}_\text{f} $

(15)

式中：n_f表示声流体耦合界面的外法线单位矢量。

此外，声固耦合面上须满足压力平衡条件

$ \sigma \left| {}_{{{\boldsymbol{n}}_\text{e}}}\right.=-p $

(16)

式中：

$ {\boldsymbol{n}}_\text{e}=-{\boldsymbol{n}}_\text{f} $

表示弹性体耦合界面的外法线单位矢量。

依据式（13）、（14）和式（16），问题域内总声压p的作用相当于对弹性体表面施加了一个外部的压力激励，同时代入散射声压的离散表达式，可得弹性力外部压力激励为

$ \begin{split}& \boldsymbol{f}_\text{e}=\left.\int_{\varGamma} \boldsymbol{\varPsi}_\text{e}^{\mathrm{T}} \boldsymbol{n}_\text{e} \sigma\right|_{\boldsymbol{n}_\text{e}} \mathrm{~d} \varGamma=\int_{\varGamma} \boldsymbol{\varPsi}_\text{e}^{\mathrm{T}} \boldsymbol{n}_\text{f} p \mathrm{~d} \varGamma=\int_{\varGamma} \boldsymbol{\varPsi}_\text{e}^{\mathrm{T}} \boldsymbol{n}_\text{f}\left(p_\text{i}+p_\text{s}\right) \mathrm{d} \varGamma =\\&\qquad \left(\int_{\varGamma} \boldsymbol{\varPsi}_\text{e}^{\mathrm{T}} \boldsymbol{n}_\text{f} \boldsymbol{\varPsi}_\text{f} \mathrm{~d} \varGamma\right) \boldsymbol{p}_\text{s}+\int_{\varGamma} \boldsymbol{\varPsi}_\text{e}^{\mathrm{T}} \boldsymbol{n}_\text{f} p_\text{i} \mathrm{~d} \varGamma=\boldsymbol{H} \boldsymbol{p}_\text{s}+\int_{\varGamma} \boldsymbol{\varPsi}_\text{e}^{\mathrm{T}} \boldsymbol{n}_\text{f} p_\text{i} \mathrm{~d} \varGamma\end{split} $

(17)

式中：H为声固耦合矩阵，

$ {\boldsymbol{\varPsi }}_\text{f} $

和

$ \boldsymbol{\varPsi }_\text{e}^{} $

分别表示流体区域和弹性体区域的形函数。

依据式（12）和式（15），弹性体表面振速（即位移导数）也会对声流体质点产生一个外部的速度激励，同时代入时谐关系式

$ {\boldsymbol{\dot{u}}}_\text{e}\left({x}_\text{i},t\right)=\text{i}\omega {\boldsymbol{\varPsi }}_\text{e}{\boldsymbol{u}}_\text{e}\left({x}_\text{i}\right) $

，可得流体质点的外部速度激励为

$ \boldsymbol{f}_\text{f}=-\mathrm{i} \rho \omega \int_{\varGamma} \boldsymbol{\varPsi}_\text{f}^{\mathrm{T}} \boldsymbol{n}_\text{f}^{\mathrm{T}} \dot{\boldsymbol{u}}_\text{f} \mathrm{~d} \varGamma=-\mathrm{i} \rho \omega \int_{\varGamma} \boldsymbol{\varPsi}_\text{f}^{\mathrm{T}} \boldsymbol{n}_\text{f}^{\mathrm{T}} \dot{\boldsymbol{u}}_\text{e} \mathrm{~d} \varGamma=\rho \omega^2\left(\int_{\varGamma} \boldsymbol{\varPsi}_\text{f}^{\mathrm{T}} \boldsymbol{n}_\text{f}^{\mathrm{T}} \boldsymbol{\varPsi}_\text{e} \mathrm{~d} \varGamma\right) \boldsymbol{u}_\text{e}=\rho \omega^2 \boldsymbol{H}^{\mathrm{T}} \boldsymbol{u}_\text{e} $

(18)

联立式（12）和式（13），并代入式（17）和式（18），可得水下弹性体声固耦合分析的控制方程

$ \left[\begin{array}{cc}\boldsymbol{K}_\text{e}-\omega^2 \boldsymbol{M}_\text{e} & -\boldsymbol{H} \\-\rho \omega^2 \boldsymbol{H}^\text{T} & \boldsymbol{K}_\text{f}-k^2 \boldsymbol{M}_\text{f}\end{array}\right]\left[\begin{array}{l}\boldsymbol{u}_\text{e} \\\boldsymbol{p}_\text{s}\end{array}\right]=\left[\begin{array}{c}\boldsymbol{F}_\text{e}+\displaystyle\int_{\varGamma} \boldsymbol{\varPsi}_\text{e}^{\mathrm{T}} \boldsymbol{n}_\text{f} p_\text{i} \mathrm{~d} \varGamma \\\boldsymbol{F}_\text{f}\end{array}\right] $

(19)

上式中，若仅考虑声辐射分析，则入射声压p_i = 0。

3　数值算例

收起

3.1　三维充液驻波管内声传播

如图3所示，考虑三维充液驻波管的声传播，其左端壁面存在一个法向速度激励v_n = 10 m/s，其余壁面设为刚性，即法向速度为0。管长为L = 1 m，宽和高分别为b = 0.1 m和h = 0.1 m。管内充水，水的密度为ρ_w = 1000 kg/m³，水中声速为c_w = 1500 m/s。通过Delaunay四面体剖分算法将问题域划分为11812个4节点四面体，作为FEM-T4的单元、以及RPIM和MRPIM的积分网格，未对网格进行额外优化，包含2640个场节点，节点平均距离d_c ≈ 0.168 m；通过Delaunay四面体剖分算法将问题域划分为1424个10节点四面体，作为FEM-T10的单元，同样未对网格进行额外优化，包含2714个场节点，节点平均距离d_c ≈ 0.166 m。图4描绘了相应的节点分布。

管内声压解析解为

$ p^{\mathrm{e}}=-\mathrm{i} \rho c v_0 \frac{\cos (k L(1-x / L))}{\sin k L} $

(20)

管内声场模态频率解析解为

$ f^{\text {eigen }}=\frac{c_{\mathrm{w}}}{2}\left(\left(\frac{m-1}{L}\right)^2+\left(\frac{n-1}{b}\right)^2+\left(\frac{o-1}{h}\right)^2\right)^{0.5} $

(21)

式中：m，n和o为与模态频率有关的非负整数参数，且不能同时为零。

首先，为保证更高的计算精度，对径向基函数中的待定参数α_c和q的取值范围进行探讨。如图5所示，α_c取0~5范围内的整数，q取0~4.0，绘制获取的数值解中的L²相对误差曲线（声压数值解p与解析解p^e之间的相对误差

${{\varepsilon }_{r}}={{\left| {{p}^{e}}-p \right|}_{{{L}^{2}}}}/{{\left| {{p}^{e}} \right|}_{{{L}^{2}}}} $

）。由图中结果可知，α_c = 0时的相对误差曲线始终最低，且q = 0.95时计算误差最小。如图6所示，q取0.95、1.03和2.83，α_c取−4.0~4.0，绘制获取的数值解中的L²相对误差曲线。由图可知，α_c = 0和q = 0.95时，数值解中的相对误差最低，因此在本文中采用了上述参数。

众所周知，内场声学问题中的声压响应可表述为各阶声学模态（与特征频率对应）的叠加。因此，精确计算声学模态是获得可靠声压数值解的关键步骤。图7展示了由不同数值方法计算得到的前80阶特征频率与式（21）给出的模态解析解之间的相对误差。相对误差定义为

$ {\varepsilon }_{m}=\left(f_i^{h}-{f}_i\right)/ {f}_i\times 100\mathrm{\% } $

，

$ f_i^{h} $

代表由特定数值方法计算得到的第i阶模态特征频率，

$ {f}_i $

代表解析解得到的第i阶模态特征频率。图中，FEM-T4表示四面体线性有限单元（4节点），FEM-T10表示四面体二阶有限单元（10节点）。

在声学FEM分析中，研究者普遍遵循一个“经验法则”^[13]，即一个波长的空间距离内至少需离散为6个单元或7个节点，相当于无量纲波数kd_c ≈ 1。遵循该法则，本算例的计算频率应不超过14800 Hz，该频率接近于第70阶模态（当m = 1, n = 1, o = 3 或 m = 1, n = 3, o = 1 时，频率为15000 Hz）。换言之，理论上采用各种数值方法应能精确计算出前70阶的模态频率。但是，使用FEM-T4得到的模态频率与解析解存在显著偏差，表明即使在经验法则适用范围内，FEM-T4数值解也难以达到合适的精度。相比之下，采用FEM-T10得到的模态频率则更接近解析解。特别的是，MRPIM和RPIM总体上能提供更精确的高阶模态频率，并且MRPIM得到的特征频率通常略优于RPIM，即相对误差更接近0%。综上所述，可以预期在此算例中MRPIM应具有最佳的计算精度。

图8展示了在1~14800 Hz（受“经验法则”约束）频率范围内，所使用的数值方法解中的L²相对误差，即声压数值解与解析解之间的相对误差。与前述预测一致，MRPIM解的数值误差普遍低于其他数值方法。值得一提的是，传统RPIM的计算精度在较低频率下下降较为明显，而在较高频率下其精度趋近于MRPIM。这可能是由于数值积分误差在低频范围内占主导，导致RPIM在较低频率下存在计算精度损失。此外，FEM-T4数值解中的相对误差对频率升高尤为敏感，且FEM-T10在14800 Hz附近时也无法保证足够的精度（其相对误差逐渐超过10%），表明这两种方法难以适应“经验法则”的要求。

本文接下来对比了MRPIM与其他数值方法在收敛性（如图9）和计算效率（如图10）方面的表现。这两项属性是评估数值方法在工程应用中可行性与实用性的基石。图9描绘了f =10000 Hz时L²相对误差随特征节点间距的变化关系，可反映不同数值方法的收敛性，C表示曲线的斜率（由于FEM-T4曲线处于“预渐进”阶段，因此未展示其斜率）。需说明的是，4节点四面体网格的节点平均距离设定为0.0769 m、0.0625 m、0.0437 m、0.0237 m、0.0168 m和0.0091 m；而10节点四面体网格的节点平均距离则设定为0.0499 m、0.0401 m、0.0227 m、0.0166 m和0.0085 m。从图中可得，随着节点平均距离的减小，MRPIM解的相对误差随之降低，逐渐收敛于解析解。此外，在相同的节点平均距离条件下，MRPIM始终表现出最低的相对误差，这意味着MRPIM只需最少的节点数即可达到预设的精度要求。

图10中描绘了相对误差与计算时间的关系，展示了在10000 Hz时不同数值方法的计算效率。计算时间（单位为秒）代表使用相同直接求解器求解系统矩阵所耗费的时间。由图可知，在所考虑的频率点处，MRPIM和FEM-T10的曲线通常位于其他方法曲线下方，这表明它们能够在更短的时间内达成更高的计算精度。综上所述，MRPIM和FEM-T10相较于FEM-T4及传统RPIM具有更优的计算效率。

3.2　水下三维弹性空心球的声散射

考虑水下三维弹性空心球在无限大水域中的声散射，弹性空心球外径r = 0.5 m，材料为钢制，球壁厚t = 0.1 m，密度ρ_e = 7850 kg/m³，杨氏模量E = 210 GPa，泊松比v = 0.3，纵波波速为c_L = 6000.9 m/s，横波波速为c_S = 3207.7 m/s。水的密度为ρ_w = 1000 kg/m³，水中声速为c_w = 1500 m/s。如图11所示，将无限大问题域缩减为6 m× 6 m× 6 m的有限域，其中最外层为厚度α = 2 m的PML吸声层，问题域中心布置了钢制空心球。通过Delaunay四面体剖分算法将问题域划分为966 704个四节点四面体，作为FEM-T4的单元和MRPIM的积分网格，未对网格进行额外优化，包含166 519个场节点。图12描绘了相应的节点分布，节点平均距离d_c ≈ 0.111 m。此问题解析解如下式

$ p(R, \theta)=\sum_{n=0}^{\infty} c_n\left[-(-\mathrm{i})^n \dfrac{J_n(k r) Q_n-k a_n \dfrac{\mathrm{~d}}{\mathrm{~d}(k r)} J(k r) S_n}{H_n^{(2)}(k r) Q_n-k a_n \dfrac{\mathrm{~d}}{\mathrm{~d}(k r)} H_n^{(2)}(k r) S_n}\right] H_n^{(2)}(k R) \cos (n \theta) $

(22)

式中：R和θ为极坐标轴方向，其余符号及参数见文献[14]。

如前所述，如遵循“经验法则”，此算例最高计算频率约为2250 Hz。在此，分析了500 Hz（kd_c ≈ 0.23）、1000 Hz（kd_c ≈ 0.47）、1500 Hz（kd_c ≈ 0.70）和2250 Hz（kd_c ≈ 1.05）下空心球的声散射，并绘制了距离球中心1.5 m处的散射声波指向性图，如图13~16所示。由图可见，FEM-T4解随着频率增长，逐渐更多地偏离了解析解，并且在1500 Hz时开始出现虚假的高频波动。这是由于FEM-T4中存在随频率快速增长的数值频散误差和数值各向异性，前者导致数值声速高于真实声速，后者导致数值声速沿传播角度呈不均匀性分布。作为对比，MRPIM解在各频率点更贴合解析解，其表现更加稳定。由此可知，需进一步使用更精细、更规则的网格以保证FEM-T4解的精度和稳定性，即声学“经验法则”不能保证FEM-T4解的有效性，而MRPIM在声学“经验法则”所规定的频率范围内表现相对良好。

为量化分析上述数值解中的误差，将不同数值方法在各频率点处的相对误差列入表1中。由表可知，FEM-T4解的相对误差在所考虑频率点处总是高于MRPIM解的相对误差。并且，在较低频率1500 Hz（kd_c ≈ 0.70）时，其相对误差已超过15%，不能满足工程中计算精度的要求，表明FEM-T4的计算精度对计算频率的高低非常敏感。反观MRPIM，其解的相对误差在2250 Hz（kd_c ≈ 1.05）时仍在10%左右，表明其在满足“经验法则”要求的频率范围内可一定程度上保证计算结果的可靠性。需要说明的是，在500 Hz（kd_c ≈ 0.23）处，MRPIM和FEM-T4解都出现了一定精度损失，这是由于所使用PML的坐标复变换与频率有关，频率越高，阻尼效果越好，因此在较低频率处的吸声效果会有所下降。

最后对MRPIM和FEM-T4的误差控制策略进行探讨，即研究每波长距离至少含几个节点可保证数值结果的可靠性。由表2可知，当每波长节点数在11~12个时，即每波长单元数在10~11个时，FEM-T4才能将声压相对误差降低至10%左右，表明工程计算中常用的“经验法则”已对FEM-T4失效。作为对比，当每波长节点数在6~7个时，MRPIM声压解的相对误差接近10%，表明“经验法则”可适用于基于MRPIM的工程计算。若按照相应的误差控制策略开展三维声学计算，在边长为一个波长的正六面体中，FEM-T4至少需要11³ = 1331个节点，MRPIM至少需要7³ = 343个节点。相对来说，MRPIM所需节点数相比FEM-T4减少了近75%，可极大降低三维声学计算中的自由度和系统矩阵规模，进而减少计算内存及存储消耗，提高计算效率。

4　结　论

收起

本文针对传统有限元法在高波数条件下因频散误差导致的精度下降问题，以及其对网格质量高度依赖的局限性，提出了一种改进径向点插值无网格计算模型（MRPIM），并将其应用于三维水下弹性目标的声散射分析。主要结论如下：

（1）MRPIM通过引入基于积分网格中心的新型插值节点选择方案，能够在单个积分单元内构造连续的声压近似函数，显著降低了数值积分误差，提高了计算精度和数值稳定性。

（2）与传统RPIM相比，MRPIM在每个积分单元内仅需进行一次“矩”矩阵求逆操作，大幅减少了计算成本，提升了计算效率。

（3）数值算例的结果表明，使用相同节点数量时，MRPIM在所考虑的中高频率范围内均表现出优于FEM的计算精度，未出现高频虚假波动现象；获取相同精度时，相比线性FEM可减少近75%节点数量，且计算效率可与二阶FEM相当，验证了其在中高频率范围内的可靠性和稳定性。

综上所述，MRPIM具有良好的工程应用前景，但目前仍缺乏在大尺度问题（如数百万级节点）中的计算成本分析、求解复杂几何模型声散射的可靠性等方面的研究等。

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

习强, 傅卓佳. 基于浅海环境潜器振动的水下声辐射计算模型研究进展[J]. 科学通报, 2022, 67(27): 3269‒3281.

, Fu

Z J

. A review of computational models for underwater acoustic radiation induced by structural vibration in the shallow marine environment[J]. China Science Bulletin, 2022, 67(27): 3269‒3281. (in Chinese)

邹明松, 吴有生. 船舶声弹性力学理论及其应用[J]. 力学进展, 2017, 47: 385‒428.

Zou

M S

, Wu

Y S

. Sano-elasticity of ships and the related applications[J]. Advances in Mechanics, 2017, 47: 385‒428. (in Chinese)

Liu

G R

, Zhang

G Y

. Smoothed point interpolation methods: G space theory and weakened weak forms[M]. Singapore: World Scientific, 2019.

Wenterodt

, Von

Estorff O

. Optimized meshfree methods for acoustics[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(25-28): 2223‒2236.

Uras

R A

, Chang

C T

, Chen

, et al. Multiresolution reproducing kernel particle methods in acoustics[J]. Journal of Computational Acoustics, 1997, 5(1): 71‒94.

Liu

G R

, Zhang

, Gu

Y T

, et al. A meshfree radial point interpolation method (RPIM) for three-dimensional solids[J]. Computational Mechanics, 2005, 36: 421‒430.

夏百战, 于德介, 姚凌云. 二维声学数值计算的径向插值有限元法[J]. 机械工程学报, 2012, 48(11): 159‒165.

Xia

B Z

, Yu

D J

, Yao

L Y

. Radial interpolation finite element method for two dimension acoustic numerical computation[J]. Journal of Mechanical Engineering, 2012, 48(11): 159‒165. (in Chinese)

姚凌云, 于德介, 臧献国. 声学数值计算的分区光滑径向点插值无网格法[J]. 振动与冲击, 2011, 30(10): 188‒192.

Yao

L Y

, Yu

D J

, Zang

X G

. Numerical computation for acoustic problems with a cell-based smoothed radial point interpolation method[J]. Journal of Vibration and Shock, 2011, 30(10): 188‒192. (in Chinese)

桂强, 孙磊, 游翔宇, 等. 二维水下声辐射问题的改进径向基点插值法研究[J]. 计算力学学报, 2024, 41(2): 344‒351.

Gui

, Sun

, You

X Y

, et al. A modified radial point interpolation method for the analysis of underwater acoustic radiation problems[J]. Chinese Journal of Computaional Mechanics, 2024, 41(2): 344‒351. (in Chinese)

吴绍维, 柯磊, 韩国文. 一种预报水下声辐射的无网格弱式径向点插值方法[J]. 振动与冲击, 2023, 42(4): 145‒155.

S W

, Ke

, Han

G W

. A meshless weak radial point interpolation method for predicting underwater acoustic radiation[J]. Journal of Vibration and Shock, 2023, 42(4): 145‒155. (in Chinese)

You

X Y

, Gui

, Zhang

Q F

, et al. Meshfree simulations of acoustic problems by a radial point interpolation method[J]. Ocean Engineering, 2020, 218(1): 108202.

You

X Y

, Chai

Y B

, Li

. Edged-based smoothed point interpolation method for acoustic radiation with perfectly matched layer[J]. Computers & Mathematics with Applications, 2020, 80(6): 1596‒1618.

Marburg

. Six boundary elements per wavelength: Is that enough?[J]. Journal of Computational Acoustics, 2002, 10(1): 25‒51.

汤渭霖, 范军, 马忠成. 水中目标声散射[M]. 北京: 科学出版社, 2018.

Tang

W L

, Fan

, Ma

Z C

. Acoustic scattering from underwater targets[M]. Beijing: Science Press, 2018. (in Chinese)

2026年第30卷第4期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.3969/j.issn.1007-7294.2026.04.013

接收时间：2025-11-11
首发时间：2026-07-07
出版时间：2026-04-15

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-11-11

基金

作者信息

^1.长沙理工大学水利与海洋工程学院，长沙 410114

^2.上海交通大学海洋工程全国重点实验室，上海 200240

^3.水沙科学与水灾害防治湖南省重点实验室，长沙 410114

^4.华中科技大学船舶与海洋工程学院，武汉 430074

通讯作者:

姚　宇（1982–），男，博士，教授，通讯作者，E-mail：yaoyu821101@163.com

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/cblx/CN/10.3969/j.issn.1007-7294.2026.04.013

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT