0 引　言

以基本解为基础的边界元方法（BEM）在求解无限或半无限域声学问题时具备天然优势^[1]，然而单独使用边界元只能计算边界上声压 $ p $ 和声通量 $ q $ 之间有确定关系的情况，如软、硬和阻抗边界。当考虑声振耦合问题时，边界上的声压和声通量不仅受到边界积分方程的制约，还必须满足结构振动方程，此时边界元法必须与其他力学方法联合求解。在过去数十年间，有限元法（FEM）作为一种有效且稳定的数值算法，在结构力学分析中得到了广泛应用。边界元法与有限元法耦合求解算法（BEM-FEM）由Everstine^[2]首次提出后，经过众多学者的发展和完善，在外场结构声学领域得到了广泛的应用^[3-5]。传统有限元法的计算精度受网格类型和质量的影响，尽管许多自适应算法已经提出，但其效率依赖于误差估计策略和网格细化算法，因此学界从20世纪90年代兴起无网格法的研究热潮^[6]。

与传统有限元方法不同，无网格法通过散布在域和边界上的一组节点来表示求解区域，通过局部支持域所包含的节点信息构造形函数。因此无网格法可以部分甚至彻底消除对网格的需求。此外，无网格法形函数的选取具备一定的灵活性，可以根据问题类型选择不同的形函数，提高计算的精度。因此无网格法可以作为有限元方法的替代与边界元方法耦合，从而避开有限元网格划分带来的问题，同时提高计算精度。

目前关于边界元耦合无网格法方面的研究主要集中在边界条件的施加上^[7-8]，对无限声场和结构分别使用边界元法和无网格法离散这一方面研究较少。径向基点插值无网格法中（RPIM）相比其他类型无网格法优势在于：①形函数具备 $ \delta $ 函数性质，便于施加本质边界条件；②只有节点距离一个变量，构造简单；③计算结果稳定，对不同的节点分布适应性较好。因此，这里采用边界元与径向基点插值无网格方法耦合（BEM-RPIM）求解无限域弹性结构声散射问题。

1 边界元耦合径向基点插值无网格法计算模型 1.1 声散射问题的边界积分方程

目前已有文献[9]对声散射问题的边界积分方程进行了推导，这里直接给出声散射问题的常规边界积分方程（CBIE）表达式：

$\begin{aligned}[b] c({{{\boldsymbol{x}}}})p({\boldsymbol{x}}) =& \int_\varGamma {G({\boldsymbol{x}},{\boldsymbol{y}})q({\boldsymbol{y}})} {\text{d}}\varGamma ({\boldsymbol{y}}) -\\ & \int_\varGamma {\frac{{\partial G({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial n({\boldsymbol{y}})}}p({\boldsymbol{y}}){\text{d}}\varGamma ({\boldsymbol{y}})} + {p_{in}}({\boldsymbol{x}}),\quad {\boldsymbol{x}} \in \varGamma\end{aligned}。$

(1)

式中： $ G({\boldsymbol{x}},{\boldsymbol{y}}) $ 表示格林函数； $ p({\boldsymbol{x}})和q({\boldsymbol{x}}) $ 分别表示边界上 $ \varGamma $ 点 $ {\boldsymbol{x}} $ 处声压和声通量 ； $ c({\boldsymbol{x}}) $ 由点 $ {\boldsymbol{x}} $ 处的几何特征决定，若点 $ {\boldsymbol{x}} $ 处光滑， $ c({\boldsymbol{x}}) = {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2} $ ； $ {p_{in}}({\boldsymbol{x}}) $ 为点 $ {\boldsymbol{x}} $ 处的入射声压值。

将式(1)对外法向 $ n({\boldsymbol{x}}) $ 求导可得超奇异边界积分方程（HBIE）：

$ \begin{split} c({\boldsymbol{x}})q({\boldsymbol{x}}) =& \int_\varGamma {\frac{{\partial G({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial n({\boldsymbol{x}})}}q({\boldsymbol{y}}){\text{d}}\varGamma ({\boldsymbol{y}})} - \\ & \int_\varGamma {\frac{{{\partial ^2}G({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial n({\boldsymbol{x}})\partial n({\boldsymbol{y}})}}p({\boldsymbol{y}}){\text{d}}\varGamma ({\boldsymbol{y}})} {\text{ + }}\frac{{\partial {p_{in}}({\boldsymbol{x}})}}{{\partial n({\boldsymbol{x}})}}, {\boldsymbol{x}} \in \varGamma。\end{split} $

(2)

对于二维声学问题，式(1)和式(2)的各阶核函数表达式如下：

$ \begin{aligned}[b] & G({\boldsymbol{x}},{\boldsymbol{y}})=\frac{i}{4}H_0^{(1)}(kr), \\ & \frac{{\partial G({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial n({\boldsymbol{y}})}} = - \frac{{ik}}{4}H_1^{(1)}(kr)\frac{{\partial r}}{{\partial n({\boldsymbol{y}})}}, \\ & \frac{{\partial G({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial n({\boldsymbol{x}})}} = - \frac{{ik}}{4}H_1^{(1)}(kr)\frac{{\partial r}}{{\partial n({\boldsymbol{x}})}}, \\ & \frac{{{\partial ^2}G({\boldsymbol{x}},{\boldsymbol{y}})}}{{\partial n({\boldsymbol{x}})\partial n({\boldsymbol{y}})}} =\frac{{ik}}{4}H_1^{(1)}(kr)n({\boldsymbol{x}})n({\boldsymbol{y}}) +\\ & \frac{{i{k^2}}}{4}H_2^{(1)}(kr)\frac{{\partial r}}{{\partial n({\boldsymbol{x}})}}\frac{{\partial r}}{{\partial n({\boldsymbol{y}})}}。\end{aligned} $

(3)

式中： $ r = \left| {{\boldsymbol{x}} - {\boldsymbol{y}}} \right| $ 表示场点与源点之间的欧式距离； $ H_n^{(1)} $ 表示 $ n $ 阶第一类汉克尔函数。

1）解的非唯一性问题

单独使用式(1)或式(2)求解外声场问题时，在某些特定频率处会出现计算结果不稳定的问题，即解的非唯一性问题^[1]，采用Burton-Miller方法^[10]解决此问题，可表示为：

$ {\text{CBIE}} + \alpha {\text{HBIE}} = 0。$

(4)

式中：CBIE和HBIE分别表示常规边界积分方程(1)和超奇异边界积分方程(2)； $ \alpha $ 为耦合系数，针对采用的基本解形式和时间项约定，耦合系数的取值^[11]为 $ \alpha = {{ - i} \mathord{\left/ {\vphantom {{ - i} k}} \right. } k} $ 。

2）边界积分方程离散

边界积分方程描述的是连续系统，为便于计算机求解，需要将其离散化为线性方程组。选择常单元配点法对式(4)进行离散可得系统方程：

$ {\boldsymbol{Hp}} = {\boldsymbol{Gq}} + {{\boldsymbol{P}}^{in}} 。$

(5)

式中：H和G是边界元系数矩阵，阶数 $ N \times N $ ； $ {{\boldsymbol{P}}^{in}} $ 为入射波对应的向量，阶数 $ N \times 1 $ ； $ p $ 和 $ q $ 为边界上的声压和声通量向量，阶数 $ N \times 1 $ 。具体表达式如下：

$ \begin{split} & {H_{mn}} = H_{mn}^1 + \alpha H_{mn}^2，\\ & {G_{mn}} = G_{mn}^1 + \alpha G_{mn}^2，\\ & P_m^{in} = {P^{in}}({{\boldsymbol{x}}_m}) + \alpha \frac{{\partial {P^{in}}({{\boldsymbol{x}}_m})}}{{\partial {\bf{n}}({{\boldsymbol{x}}_m})}}，\end{split} $

(6)

$ \begin{split} & H_{mn}^1 = \int_{{\varGamma _n}} {\frac{{\partial G({{\boldsymbol{x}}_m},{\boldsymbol{y}})}}{{\partial n({\boldsymbol{y}})}}{\text{d}}{\varGamma _n}({\boldsymbol{y}})} + {\delta _{mn}}c({{\boldsymbol{x}}_m}) ，\\ & H_{mn}^2 = \int_{{\varGamma _n}} {\frac{{{\partial ^2}G({{\boldsymbol{x}}_m},{\boldsymbol{y}})}}{{\partial n({{\boldsymbol{x}}_m})\partial n({\boldsymbol{y}})}}{\text{d}}{\varGamma _n}({\boldsymbol{y}})} ，\\ & G_{mn}^1 = \int_{{\varGamma _n}} {G({{\boldsymbol{x}}_m},{\boldsymbol{y}}){\text{d}}{\varGamma _n}} ({\boldsymbol{y}}) ，\\ & G_{mn}^2 = \int_{{\varGamma _n}} {\frac{{\partial G({{\boldsymbol{x}}_m},{\boldsymbol{y}})}}{{\partial n({{\boldsymbol{x}}_m})}}{\text{d}}{\varGamma _n}({\boldsymbol{y}})} - {\delta _{mn}}c({{\boldsymbol{x}}_m})。\end{split} $

(7)

1.2 径向基点插值无网格法（RPIM）

采用基于径向基点插值形函数的弱形式无网格法^[12]对结构域进行离散，弱形式无网格法（见图1）与传统有限元的主要区别在于形函数的构造方式。

1）径向基点插值形函数

添加多项式的径向基点插值可表示为：

$ {\boldsymbol{u}}({\boldsymbol{x}}) = \sum\limits_{n = 1}^N {{a_n}{R_n}({\boldsymbol{x}}{\text{)}}} + \sum\limits_{m = 1}^M {{b_m}{P_m}({\boldsymbol{x}}) = {{R}^{\rm{T}}}({\boldsymbol{x}}) \cdot {\boldsymbol{a}} + {{P}^{\rm{T}}}({\boldsymbol{x}}) \cdot {\boldsymbol{b}}}。$

(8)

式中： $ {R_n}({\boldsymbol{x}}{\text{)}} $ 为径向基函数（RBF）； $ N $ 为局部主持域中包含的节点个数； $ {P}({\boldsymbol{x}}) = {[1{\text{ }}x{\text{ }}y]^{\rm{T}}} $ 表示二维线性多项式函数； $ {a_n} $ 和 $ {b_m} $ 为待定系数。

在径向基函数 $ {R_n}({x}{\text{)}} $ 中，仅有表示计算点 $ {\boldsymbol{x}} $ 和节点 $ {{\boldsymbol{x}}_{{n}}} $ 之间距离的变量为：

$ r = \left| {{\boldsymbol{x}} - {{\boldsymbol{x}}_n}} \right|，$

(9)

采用复合2次（MQ）函数，即

$ {R_n}(r) = {[{r^2} + {({a_c}{d_c})^2}]^q} 。$

(10)

其中： $ {a_c} $ 和 $ q $ 是MQ-RBF的2个形状参数，可由数值试验确定；当 $ q = \pm 0.5 $ 时为标准MQ-RBF，对于二维固体问题， $ q $ 取0.98或1.03可取得良好的分析结果^[12]； $ {d_c} $ 为平均节点间距。

为确定式(8)中的待定系数 $ {a_n} $ 和 $ {b_m} $ ，需形成计算点 $ {x} $ 的支持域，其中包含 $ N $ 个场节点，代入式(8)得到

$ {\boldsymbol{U}} = {{\boldsymbol{R}}_0}{\boldsymbol{a}} + {{\boldsymbol{P}}_m}{\boldsymbol{b}}，$

(11)

式中， $ {\boldsymbol{U}} $ 为节点位移向量。

$ {\boldsymbol{U}} = {[{u_1}{\text{ }}{u_2}{\text{ }} \cdots {\text{ }}{u_N}]^{\rm{T}}} 。$

(12)

RBFs的矩量矩阵为：

$ {{\boldsymbol{R}}_0} = {\left[ {\begin{array}{*{20}{c}} {{R_1}\left( {{r_1}} \right)}&{{R_2}({r_1})}& \cdots &{{R_N}({r_1})} \\ {{R_1}({r_2})}&{{R_2}({r_2})}& \cdots &{{R_N}({r_2})} \\ \vdots & \vdots & \ddots & \vdots \\ {{R_1}({r_N})}&{{R_2}({r_N})}& \cdots &{{R_N}({r_N})} \end{array}} \right]_{\left( {N \times N} \right)}}{\text{ }} 。$

(13)

多项式矩量矩阵为：

$ {{\boldsymbol{P}}_m} = {\left[ {\begin{array}{*{20}{c}} 1&{{x_1}}&{{y_1}} \\ 1&{{x_2}}&{{y_2}} \\ \vdots & \vdots & \vdots \\ 1&{{x_N}}&{{y_N}} \end{array}} \right]_{N \times 3}} 。$

(14)

RBFs的待求系数向量为：

$ {\boldsymbol{a}} = {[{a_1}{\text{ }}{a_2}{\text{ }} \cdots {\text{ }}{a_N}]^{\rm{T}}} 。$

(15)

多项式系数向量为：

$ {\boldsymbol{b}} = {[{b_1}{\text{ }}{b_2}{\text{ }} \cdots {\text{ }}{b_M}]^{\rm{T}}}。$

(16)

式(13)中， $ {R_i}({r_k}) $ 的 $ {r_k} $ 定义为：

$ {r_k} = \sqrt {{{({x_k} - {x_i})}^2} + {{({y_k} - {y_i})}^2}} 。$

(17)

式(11)中有 $ M + N $ 个变量，可使用下面约束条件再添加 $ M $ 个方程。

$ \sum\limits_{n = 1}^N {{P_m}({x_n}){a_n} = {\boldsymbol{P}}_m^T{\boldsymbol{a}}} = {\boldsymbol{0}},\quad m = 1,2, \cdots ,M 。$

(18)

联立式(11)和式(18)可得到方程：

$ \widehat {\boldsymbol{U}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{U}} \\ {\boldsymbol{0}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{R}}_0}}&{{{\boldsymbol{P}}_m}} \\ {{\boldsymbol{P}}_m^{\rm{T}}}&{\boldsymbol{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\boldsymbol{a}} \\ {\boldsymbol{b}} \end{array}} \right] = {\boldsymbol{G}}\widehat {\boldsymbol{a}} 。$

(19)

因为矩阵 $ {{\boldsymbol{R}}_0} $ 是对称的，故矩阵 $ {\boldsymbol{G}} $ 也是对称的。求解上式可得：

$ \widehat {\boldsymbol{a}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{a}} \\ {\boldsymbol{b}} \end{array}} \right] = {{\boldsymbol{G}}^{ - 1}} {\boldsymbol{U}} 。$

(20)

将(20)代入(8)可得：

$ {\boldsymbol{u}}({\boldsymbol{x}}) = {{R}^{\rm{T}}}({\boldsymbol{x}}) \cdot {\boldsymbol{a}} + {{P}^{\rm{T}}}({\boldsymbol{x}}) \cdot {\boldsymbol{b}} = \left[ {{{R}^{\rm{T}}}({\boldsymbol{x}}){\text{ }}{{P}^{\rm{T}}}({\boldsymbol{x}})} \right]{{\boldsymbol{G}}^{ - 1}} {\boldsymbol{U}} = { {\boldsymbol{\varPhi }}^{\rm{T}}}({\boldsymbol{x}}){\boldsymbol{U}} 。$

(21)

式中，RPIM形函数为

$ \begin{split} {{\boldsymbol{\varPhi }}^T}({\boldsymbol{x}}) & = \left[ {{{R}^{\rm{T}}}({\boldsymbol{x}}){\text{ }}{{P}^T}({\boldsymbol{x}})} \right]{{\boldsymbol{G}}^{ - 1}} =\\ & \left[ {{\phi _1}({\boldsymbol{x}}){\text{ }}{\phi _2}({\boldsymbol{x}}){\text{ }} \cdots {\text{ }}{\phi _N}({\boldsymbol{x}}){\text{ }}{\phi _{N + 1}}({\boldsymbol{x}}){\text{ }} \cdots {\text{ }}{\phi _{N + M}}({\boldsymbol{x}})} \right] 。\end{split} $

(22)

最后RPIM形函数表示为：

$ {\boldsymbol{\varPhi }}({\boldsymbol{x}}) = {\left[ {{\phi _1}({\boldsymbol{x}}){\text{ }}{\phi _2}({\boldsymbol{x}}){\text{ }} \cdots {\text{ }}{\phi _N}({\boldsymbol{x}})} \right]^{\rm{T}}} ，$

(23)

式(21)可写为：

$ {\boldsymbol{u}}({\boldsymbol{x}}) = {{\boldsymbol{\varPhi }}^{\rm{T}}}({\boldsymbol{x}}){\boldsymbol{U}} ，$

(24)

对(24)求偏导可得：

$ {{\boldsymbol{u}}_{,l}}({\boldsymbol{x}}) = {\boldsymbol{\varPhi }}_{,l}^{\rm{T}}({\boldsymbol{x}}){\boldsymbol{U}} 。$

(25)

式中：逗号表示对其后的空间坐标系求偏导数； $ l $ 表示坐标 $ x $ 或 $ y $ （二维问题）。

2）无网格法弹性结构方程

为方便积分运算，以有限元网格作为背景积分网格，通过上述方式可以构造积分点处RPIM形函数及其偏导数，之后按照有限元法分析流程得到弹性结构简谐振动方程：

$ ( - {\omega ^2}{\boldsymbol{M}} + {\boldsymbol{K}}){\boldsymbol{u}} = {\boldsymbol{F}} ，$

(26)

$ \begin{split}&{K}_{e}={\displaystyle {\int }_{{{S}}_{e}}{B}_{{{S}}_{e}}^{\text{T}}D{B}_{{{S}}_{e}}{{{\rm{d}}S}}_{e}}，\\ &{M}_{e}={\displaystyle {\int }_{{{S}}_{e}}{\rho }_{S}{N}_{{{S}}_{e}}^{\text{T}}{N}_{{{S}}_{e}}{{{\rm{d}}S}}_{e}}，\\ &K=\underset{{e=}1}{\overset{{M}}{\cup }}{K}_{e}\text{，}M=\underset{{e=}1}{\overset{{M}}{\cup }}{M}_{e}。\end{split} $

(27)

式中： $ \omega $ 为角频率， $ {\rho _S} $ 为结构密度， $ {\boldsymbol{u}} $ 为位移向量， $ {\boldsymbol{M}} $ 和 $ {\boldsymbol{K}} $ 分别为结构的整体质量矩阵和刚度矩阵， $ {\boldsymbol{F}} $ 为节点载荷矩阵； $ {{\boldsymbol{N}}_{{{\text{S}}_e}}} $ 为插值函数矩阵， $ {{\boldsymbol{B}}_{{{\text{S}}_e}}}{\text{ = }}\partial {{\boldsymbol{N}}_{{{\text{S}}_e}}} $ 为应变矩阵， $ {\boldsymbol{D}} $ 为弹性系数矩阵， $ \cup $ 表示单元的组装过程。

1.3 边界元和径向基点插值无网格法耦合方程

对于无限域弹性结构声散射问题，边界声压向量 $ p $ 作用于结构表面，在没有其他激励的情况下

$ {\boldsymbol{F}} = {\boldsymbol{Lp}} ，$

(28)

$ {{\boldsymbol{L}}_e} = \int_{{\varGamma _e}} {{{\boldsymbol{N}}_{{S_e}}}{{\boldsymbol{n}}_{{\varGamma _e}}}{\boldsymbol{N}}_{_{{F_e}}}^{\text{T}}{\text{d}}{\varGamma _e}}，$

(29)

$ {\boldsymbol{L}} = \mathop \cup \limits_{{{e = }}1}^{{M}} {{\boldsymbol{L}}_e}。$

(30)

式中： $ {\boldsymbol{L}} $ 为结构流场之间的耦合矩阵； ${{\boldsymbol{N}}_{{{{F}}_e}}}$ 和 ${{\boldsymbol{N}}_{{{{S}}_e}}}$ 分别为声场和结构的插值形函数的矩阵形式； $ {{\boldsymbol{n}}_{{\Gamma _e}}} $ 为边界单元的单位法向量。

有限元位移 $ {\boldsymbol{u}} $ 与声通量 $ {\boldsymbol{q}} $ 之间的关系为：

$ {\boldsymbol{q}} = {\omega ^2}{\rho _f}{{\boldsymbol{\Theta }}^{ - 1}}{{\boldsymbol{L}}^{\text{T}}}{\boldsymbol{u}}，$

(31)

$ {{\boldsymbol{\Theta }}_e} = \int_{{\varGamma _e}} {{{\boldsymbol{N}}_{{{{F}}_e}}}{\boldsymbol{N}}_{{{{F}}_e}}^{\text{T}}{\text{d}}{\varGamma _e}} ，$

(32)

$ {\boldsymbol{\Theta }} = \mathop \cup \limits_{{{e = }}1}^{{M}} {{\boldsymbol{\Theta }}_e} 。$

(33)

式中， $ {\boldsymbol{\Theta }} $ 为边界质量矩阵。

将式(28)和式(31)分别代入无网格法系统方程(26)和边界元方程(5)并联立得耦合系统方程：

$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}{{ - }}{\omega ^2}{\boldsymbol{M}}}&{ - {\boldsymbol{L}}} \\ {{{ - }}{\omega ^2}{\rho _f}{\boldsymbol{G}}{{\boldsymbol{\Theta }}^{ - 1}}{{\boldsymbol{L}}^{\rm{T}}}}&{\boldsymbol{H}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\boldsymbol{u}} \\ {\boldsymbol{p}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{0}} \\ {{P^{in}}} \end{array}} \right] 。$

(34)

求解上述方程组可以得到边界声压 $ p $ 和结构节点位移 $ {\boldsymbol{u}} $ ，通过式(2.31)得到边界声通量 $ {\boldsymbol{q}} $ ，之后通过域内积分方程，可求得任意一点的散射声压 $ {{\boldsymbol{p}}_s} $ ，即

$ {{\boldsymbol{p}}_s} = - {{\boldsymbol{H}}_f}{\boldsymbol{p}} + {{\boldsymbol{G}}_f}{\boldsymbol{q}} 。$

(35)

2 数值算例和结果讨论

以无限域中平面波入射无限长弹性圆柱体散射问题（可简化为平面问题）为例，通过结构模态频率计算，验证径向基点插值无网格法（RPIM）相比传统有限元法（FEM）计算精度高的特点。从频域散射声压计算精度、收敛性^[13]和对不规则网格的适应性等3个方面比较边界元耦合径向基点插值无网格法（BEM-PRIM）和传统边界元耦合有限元法（以下简称BEM-FEM）的性能。

2.1 计算模型介绍

如图2所示，平面波垂直入射浸没在无限大水域中的无限长弹性圆柱体，圆柱体半径 $r = 1\;{\text{m}}$ ，计算距离其几何中心 $R= 1\;000\;{\text{m}}$ 处，沿周向均匀分布的场点 $ xp $ 处的散射声压 $ {p_s} $ ，具体参数设置见表1。

入射声压为 $ 1\;{\text{Pa}} $ 的平面波时，该模型声散射解析解为^[14]：

$ {p}_{s}(R,\varphi )={\displaystyle \sum _{n=0}^{\infty }{\epsilon }_{n}{(1i)}^{n}{b}_{n}{H}_{n}^{(1）}(kR)\mathrm{cos}(n\phi )}。$

(36)

其中： $ {\varepsilon _n} $ 为Neumann因子； $ k $ 为入射波波数； $ {H}_{n}^{(1）} $ 为第一类 $ n $ 阶汉克尔函数； $ {b_n} $ 为散射系数，由结构振动边界条件确定；无穷级数最高项数取 $ {n_{\max }} \geqslant ka + 5 $ 即可保证计算精度^[15]，其中 $ a $ 是计算目标的特征尺寸，该问题中 $ a = 1\;{\text{m}} $ 。

采用相对误差 $ RE $ 和平均相对误差 $ ARE $ ^[15]描述数值方法的计算精度，定义如下：

$ R{E_m} = \frac{{\left| {{{\bar p}_m} - {p_m}} \right|}}{{\left| {{{\bar p}_m}} \right|}},\quad m = 1,2, \cdots ,NT ，$

(37)

$ ARE = \sqrt {{{\sum\limits_{m = 1}^{NT} {{{\left| {{p_m} - {{\bar p}_m}} \right|}^2}} } \mathord{\left/ {\vphantom {{\sum\limits_{m = 1}^{NT} {{{\left| {{p_m} - {{\bar p}_m}} \right|}^2}} } {\sum\limits_{m = 1}^{NT} {{{\left| {{{\bar p}_m}} \right|}^2}} }}} \right. } {\sum\limits_{m = 1}^{NT} {{{\left| {{{\bar p}_m}} \right|}^2}} }}} 。$

(38)

式中： $ {\bar p_m} $ 表示计算点处的解析值或高精度的数值计算值（通过细化网格得到）； $ {p_m} $ 表示数值计算值； $ NT $ 表示计算场点的数量。

2.2 RPIM与FEM结构模态频率计算误差比较

以无限长弹性圆柱结构为研究对象，对RPIM和FEM在计算模态频率精度方面进行比较。前20阶模态频率计算结果和相对误差如表2所示，模态频率误差分布如图3所示。

可以看出：

1） 计算结构模态频率时，RPIM无网格法相比FEM计算精度更高；

2） 随着模态阶数增加，FEM计算误差有增加的趋势，而RPIM计算误差仍然维持在较低水平；

3） FEM计算模态频率相比参考值均偏大，体现出“过硬”的特点；而RPIM计算模态频率围绕参考值上下波动，说明RPIM对系统刚度有适当的“软化”效果。

2.3 频域散射声压计算误差分析

按照有限元计算声场问题时，每个波长范围内划分不少于6个单元的经验准则^[16]，使用的网格尺寸可以计算的无量纲频率取值范围是 $ ka \leqslant 12 $ ，该例取 $ ka \in [1,16] $ 。场点 $ xp $ 处的散射声压 $ {p_s} $ 平均相对误差 $ ARE $ 随波数变化曲线如图4所示。

可以看出：

1） BEM-FEM和BEM-RIPM的计算误差均随着计算波数增加而增加，这是因为在计算不同波数时采用的是同一套节点和网格数据，随着波数增加，每个波长范围内划分的单元数量相应减少，因此计算误差增加；

2） 在波数较低的情形下（ $ ka < 6 $ ），2种方法计算精度相当且较高（ $ {\text{ARE < 1}}{{\text{0}}^{ - 2}} $ ）；但是在此之后（ $ ka > 6 $ ），BEM-RPIM计算精度明显优于BEM-FEM。

2.4 收敛性分析

为了比较BEM-RPIM和BEM-FEM的收敛性，以无量纲频率 $ka{\text{ = 3}}\text{π} {\text{,5}}\text{π}$ 为例，通过逐步细化网格，节点自由度（DOF）从80增加到7428，分别计算2种方法在场点处散射声压的平均相对误差，如图5所示。

可以看出：

1） 随着网格的细化，2种方法在计算频率处的平均相对误差均下降，且下降速率基本一致，这说明BEM-RPIM具备和BEM-FEM相似的收敛速度；

2） 在计算同一频率（ $ka = 3\text{π}$ 或 $5\text{π}$ ），自由度相同时，BEM-RPIM平均相对误差更小，即在保证相同计算精度时，BEM-RPIM需要更少的节点自由度，需要的计算资源更少。

2.5 不规则网格适应性分析

无网格方法相比有限元方法的一大优势是其形函数构造不依赖预先定义的网格，因此BEM-RPIM只需对问题域边界划分网格，相比对整个问题域划分网格工作量明显下降。而且RPIM形函数由散布在支持域内的节点定义，相比传统FEM对网格质量要求不高。为了说明BEM-RPIM同样具有这一特点，采用规则网格和不规则网格对比散射声压计算精度。按照文献[17]采取在原有节点施加随机扰动的方式获得不规则网格。

$ {\boldsymbol{x'}} = {\boldsymbol{x}} + {d_c} \cdot {r_c} \cdot \beta 。$

(39)

式中： $ {\boldsymbol{x'}} $ 为不规则网格的节点坐标； $ {\boldsymbol{x}} $ 为规则网格的节点坐标； $ {d_c} $ 为平均节点间距； $ {r_c} $ 为随机变化系数， $ {r_c} \in [ - 1,1] $ ； $ \ \beta $ 表征不规则网格的扭曲程度， $ \ \beta \in [0,0.5] $ 。 $ \ \beta $ 取值越大表示网格越扭曲。以本文所取的 $ \ \beta = 0.4 $ 为例，图6展示了规则和不规则网格的对比情况。

选取无量纲波数 $ ka = 3{\text{π}} ,5{\text{π}} $ ，按分布角度 $ \phi \in [0,{180^\circ }] $ 绘制场点 $ xp $ 处的散射声压相对误差变化如图7和图8所示。

可以看出：

1） 散射声压相对误差沿场点分布角度有比较大的波动，这是因为在某些角度散射声压参考解 $ \left| {{{\bar p}_s}} \right| $ 较小，数值解很小的偏差（ $ \left| {{{\bar p}_s} - {p_s}} \right| $ ）就会使得相对误差较大。

2） 不规则网格时的计算结果相比规则网格整体相对误差变大，这说明不规则网格对2种方法计算精度均会产生不利影响。虽然无网格方法形函数构造不依赖网格，但依赖支持域内的节点，构造RPIM形函数时采用和有限元同一套节点信息，生成不规则网格时改变了节点的分布，一定程度上破坏了节点分布的均匀性，因此导致BEM-RPIM计算精度同样下降。但值得注意的是：采用不规则网格的BEM-RPIM方法在大部分场点的计算精度仍优于采用规则网格的BEM-FEM。

3 结　语

本文将边界元方法和径向基点插值无网格方法结合用以求解无限域弹性结构声散射问题。综合对BEM-RPIM和BEM-FEM在计算精度、收敛性和对不规则网格适应性的分析比较，可以得出如下结论：

1） RPIM无网格法相比传统有限元在计算结构模态频率方面具备更高的计算精度。

2） BEM-RPIM具备和BEM-FEM相似的收敛速率，但BEM-PRIM数值精度更高：在相同节点自由度下BEM-RPIM计算精度更高；在保证相同精度时，BEM-RPIM相比BEM-FEM需要的节点自由度更少。

3） BEM-RPIM相比BEM-FEM对网格质量的依赖性更弱，虽然不规则的节点或网格分布同样会降低BEM-RPIM的计算精度，但是BEM-RPIM在不规则网格下的表现依然优于BEM-FEM在规则网格下的表现。