预报振动噪声的径向基点插值无网格与无限元耦合方法_吴绍维.pdf

 振动与冲击 第 39 卷第 10 期JOURNAL OF VIBRATION AND SHOCKVol．39 No．10 2020 基金项目国家自然科学基金项目 51909016 ; 重庆市教育委员会科学 技术研究项目 KJQN201900705 收稿日期2018 －10 －12修改稿收到日期2019 －01 －28 第一作者 吴绍维 男， 博士， 讲师， 1987 年生 通信作者 向阳 女， 博士， 教授， 1962 年生 预报振动噪声的径向基点插值无网格与无限元耦合方法 吴绍维1，向阳2，黄庭瑞3 1． 重庆交通大学航运与船舶工程学院， 重庆 400074; 2． 武汉理工大学能源与动力工程学院， 武汉 400063; 3． 安庆中船柴油机有限公司， 安徽安庆246005 摘要为克服传统的有限元耦合无限元方法中的单元匹配问题， 研究了径向基点插值法和无限元法的耦合规 律， 提出了一种预报无限域结构振动噪声的径向基点插值无网格与可变阶无限声波包络单元耦合方法， 推导了预报声压 的计算公式。为提高声场预报精度和满足声波在无限域的自由衰减， 结构外部无限声场分为使用无网格表示的近场和可 变阶声波包络单元离散的远场。在该耦合方法中， 通过在近场与远场之间的交界面上配置虚拟网格来构造具有连续性的 声压形函数， 确保了声压的连续与一致性。采用数值仿真和试验对该耦合方法进行了验证， 结果表明该耦合方法拥有无 网格法的高精度和可变阶声波包络单元法满足声波自由衰减的特点， 具有良好的精度和收敛性， 可用于实际噪声预报。 关键词噪声; 径向基点插值法 RPIM ; 声波包络单元 WEE ; 无限域 中图分类号O422． 6; O429文献标志码ADOI 10． 13465/j． cnki． jvs． 2020． 10． 005 A meshless radial point interpolation coupled with improved infinite element for predicting sound radiation WU Shaowei1，XIANG Yang2，HUANG Tingrui3 1． School of Shipping and Naval Architecture，Chongqing Jiaotong University，Chongqing 400074，China; 2． School of Energy and Power Engineering，Wuhan University of Technology，Wuhan 400063，China; 3． Anqing CSSC Diesel Engine Co． ，Ltd． ，Anqing 246005，China Abstract To eliminate the element matching in traditional based on finite and infinite elements，the coupling of the radial point interpolation RPIMand the variable order infinite acoustic wave envelope element WEE was studied． A RPIM coupled with improved WEE was proposed for exterior sound problems in infinite domain． The coupling ula of this was theoretically derived． To improve accuracy and properly imitate the amplitude decay of the outgoing travelling wave，the infinite acoustic field was divided into near and far acoustic fields， which were discretized by arbitrarily distributed nodes and WEEs respectively． The continuity and compatibility of the acoustic pressure were maintained by constructing hybrid acoustic pressure shape function on the fictitious meshin the near field along the interface between the two fields． Simulations and experiments were pered to test the ． The results show that the can take full advantage of both the RPIM and WEE s and is a valuable approach． Key wordssound;radial point interpolation RPIM ;infinite acoustic wave envelope element WEE ; infinite domain 在预报无限域结构振动噪声领域， 边界元法 Boundary Element ，BEM 是最为常用的方法， 该方法满足无穷远处的 Sommerfeld 辐射条件， 数值求 解只需离散结构边界表面， 相比最初的 d 维问题域， 目 标声场量的求解是在 d- 1 维问题域上进行， 从而降低了 问题域的维数并形成相对小的系统方程［1 －2 ］。但是在 常规的 BEM 中， 形成的系统矩阵是非稀疏和非对称矩 阵， 这显著增加了数据处理时间和储存要求， 抵消了d- 1 维问题域带来的计算效率优势。除此之外， 在与内部 问题对应的特征频率处还将出现非唯一性问题。为克 服 BEM 中的非唯一性问题， Combined Helmholtz Inte- gral Equation ulation CHIEF 方法和 Burton- Miller 方法是最常用的两种方法［3 ］。 CHIEF 点法是在结构边界节点沿内法向处添加一 ChaoXing 定数量的点 称为 CHIEF 点 ， 在这些 CHIEF 点处形成 Helmholtz 积分方程， 并与原边界 Helmholtz 积分方程联 立形成超定方程， 最后求解出最小二乘意义下的唯一 解。在这种方法中， 选取的 CHIEF 点决定该方法的稳 定性， 当选取的 CHIEF 点与结构边界面为界的内部 Dirichlet 问题的振型节点重合时， 这种方法失效。关于 如何布置 CHIEF 点和确定适用的频段， 目前还没有严 格的准则。Burton- Miller 法是在场点处对 Helmholtz 方 程求导得到超奇异积分方程， 并与原 Helmholtz 方程联 立来求解声场。使用这种方法可在全频段有效解决非 唯一性问题， 但是引入了格林函数各类奇异积分问题， 这些奇异积分需要进行复杂的数学处理， 导致计算效 率显著下降。在过去的几十年里， BEM 的巨大发展和 进步正推动着 BEM 在无限域声辐射预报中的应用， 如 近几年发展起来的快速多极 BEM 显著提高了计算效 率。目前， BEM 的相关工作仍在研究和进行中［4 －5 ］。 使用有限元法 Finite Element ，FEM 预报 无限域振动噪声面临两个问题 ［6 ］ 处理无限域并满足 无穷远处的 Sommerfeld 辐射条件; 计算频率升高引起 的数值色散问题。与边界元法不同， 使用有限元法时， 需采用人工边界将无限域截断为有限计算域， 并在人 工边界处施加无反射边界条件来代替无穷远处的 Som- merfeld 辐射条件。由于无限域被截断为有限域， 计算 时间大为减少。对于处理无限域的声辐射问题， 吸收 边界条件、 Dirichlet to Neumann DtN边界条件和完美 匹配层已被证明是有效的无反射边界条件。在这些无 反射边界条件中， DtN 边界条件建立了 Dirichlet 量与 Neumann 量之间的积分解析关系， 是一种理论上精确 的无反射边界条件［7 ］。需要注意的是 有限域越大， 消 除反射声波的效果越好， 但计算耗时增加， 因此人工边 界的选取需平衡声场预报精度和计算效率。 在基于网格的数值方法中， 随着计算频率的升高， 数值解的精度随之下降， 当计算频率超过一定阈值， 数 值解的精度将急剧下降。一种直接提高数值解精度的 措施是对网格进行细化， 但是计算时间急剧增加。为 抑制这种数值色散误差， 一些基于光滑技术或无网格 方法软化有限元模型中刚度矩阵 “过硬” 性质的方法被 提出， 如在新发展的有限元方法领域中， 为求解固体力 学问题中的偏微分方程， Liu 等 ［8 ］基于 G 空间理论提出 了新的弱化形式来解决传统有限元法中刚度矩阵“过 硬” 的问题。基于 W2 形式和源于无网格的梯度光滑 技术， 一些求解固体力学和声学的优秀方法被提出， 如 基于无网格单元的平滑 α 径向点插值方法 ［9 ］、 混合光 滑有限元法 ［10 ］、 基于稳定节点的平滑有限元方法［11 ］、 基于边缘的平滑有限元方法［12 ］以及超收敛 α 有限元 方法 ［13 ］。相比系统矩阵 “过硬” 的 FEM， 使用这些方法 和技术对系统刚度矩阵进行软化， 使其接近真实的系 统矩阵， 能够获得更精确的数值解。 20 世纪 90 年代提出的可变阶无限声波包络单元 Variable order Infinite Acoustic Wave Envelope Ele- ment，WEE 法是一种有效模拟声波幅值在无限域传 播衰减的数值技术 ［14 －16 ］， 其优点是 能够很简易地处 理无限域， 不存在格林函数各类奇异积分和声场预报 的非唯一性问题。基于并行计算策略， 该方法的计算 效率得以显著提升。但是作为一种数值方法， WEE 方 法也存在应用的局限性， 对于几何外形复杂的结构， WEE 需 与 FEM 耦 合 来 计 算 结 构 向 外 辐 射 的 声 场 ［17 －18 ］。在这种耦合方法中， WEE 与有限元单元是通 过单元匹配进行连接。当对有限元单元 或 WEE 进 行细化或改变以提高声场预报精度时， 需对 WEE 或 有限元单元 进行重构并且需再次组装系统方程， 这将 增加前处理和计算时间。由于采用了有限元单元， 该 耦合方法还存在数值色散问题。 无网格法在数值计算领域是一种比基于网格的方 法相对新的数值方法。根据公式导出的方法可将无网 格法分为三类 ［19 ］ 基于弱式的无网格法 如 Element- Free Galerkin ［20 ］ 、 基于配点技术的无网格法 如 Meshfree Collocation with Radial Basis Func- tion［21 ］ 和基于弱式和配点技术相结合的无网格法［22 ］。 无网格法是利用一组散布在问题域和域边界上的场节 点表示 而非离散 该问题域及其边界， 并利用这些任 意分布的场节点构造场变量形函数， 其本质特点是不 需要事先定义的网格或节点连接信息用于构造场变量 未知函数的插值或近似表达式。在无网格法中， 移动 最小二乘近似、 单位分解法以及径向基点插值法常用 来构造形函数。 点插 值 法 Meshfree Point Interpolation ， PIM 是 Liu 等 ［23 ］基于加辽金弱式提出的一种无网格 插值技术， 该方法基于局部分布的节点来构造形函数， 并采用一组背景网格来计算加辽金弱式积分。PIM 的 主要特点是其形函数拥有 Kronecker delta 函数性质， 可 以像 FEM 一样很容易地施加本质边界条件。现有两 类 PIM， 分别为多项式基型以及径向基型， 即多项式点 插值法以及径向基点插值法。在多项式点插值法中， 矩矩阵可能是奇异的。使用两阶矩阵三角化算法可克 服该问题 ［24 ］， 该算法能自动排除用于形成矩矩阵的节 点和多项式基的项。但由于多项式 PIM 形函数的不相 容性， 基于 Galerkin 弱式的多项式 PIM 对于不规则节 点分布不具有鲁棒性。另外， 由于局部支持域中使用 过多节点而产生的高阶多项式会导致多项式 PIM 形状 函数的变化过于剧烈。 径向基点插值法 Radial Point Interpolation Meth- 33第 10 期吴绍维等预报振动噪声的径向基点插值无网格与无限元耦合方法 ChaoXing od，RPIM 采用径向基函数 ［25 ］作为基函数构造形函 数 ［26 ］， 其对任意的节点分布具有良好的稳定性和鲁棒 性， 能够克服多项式 PIM 中的奇异性问题， 并具有比有 限元法更快的收敛率。该方法已被成功应用于解决诸 多问题， 如二维和三维固体力学、 土木工程中的材料非 线性问题 ［27 ］、 智能材料问题以及板壳结构等问题［28 ］。 关于 Helmholtz 方程的求解， RPIM 在求解内部声学问 题时已展现出高精度的优势［29 ］。若直接将 RPIM 应用 于无限域中的外部声辐射问题， 需使用无穷节点来表 示外部无限大声场， 这使得前处理以及后续计算无法 实现， 无网格法在无限域声学问题中还未得到广泛 应用 ［30 ］。 鉴于 WEE 能有效模拟声波在无限域中的传播衰 减， 将 RPIM 与 WEE 进行耦合以利用各自的优点同时 避免其缺点来计算无限域中的声学问题值得探索和研 究。本文研究了 RPIM 与 WEE 的耦合规律， 基于改进 型 WEE， 提出了预报无限域结构振动噪声的径向基点 插值无网格与无限元耦合方法， 详细给出了该方法的 公式推导和实施。在该方法中， 结构外部的无限声场 分为采用一组场节点表示的近场和 WEE 离散的远场， 通过在近场与远场之间的交界面上配置虚拟网格， 构 造具有连续性的声压形函数来完成 RPIM 与 WEE 之间 的耦合和满足声压的连续与一致性。最后采用仿真和 试验对所提出的方法进行了验证。 1理论背景 在无限大自由场 V 中， 结构做稳态简谐振动辐射 的声压满足如下形式的 Helmholtz 方程 ［31 ］ Δ 2p x k2p x 0，x ∈ V 1 式中k ω/c 为波数;ω 为辐射体边界表面振动角频 率;c 为声音在流体介质中的传播速度。声压在振动 结构边界表面 Sb上满足 Δ p x nb － ikρcv x 2 式中i2 －1;v x 为结构边界法向速度;nb 为结构 边界法向。声压应满足的 Sommerfeld 辐射条件为 lim r→∞ rα p r r ikp r 0 3 式中r 为场点到坐标原点的距离;α 为维度参数， 对 二维和三维问题分别取 1 和 2。在无穷远处， 此辐射条 件与 ρc 阻抗条件等效， 即 Δ pn∞ ikp 0， x∈S∞。 采用一虚拟交界面 Sγ将结构外部的无限大区域 V 分为内区域 Vinner和外区域 Vouter， 如图 1 所示。内区域 Vinner采用一组任意分布点场节点表示， 外区域 Vouter采 用一层 WEE 离散。假设声压在 Sγ上满足 p n γ － ikρc p Z 4 图 1结构外部无限域划分示意图 Fig． 1Division of infinite domain surrounding radiator 式中nγ为外区域 Vouter在 Sγ上的单位外法向， 其方 向指向结构;Z 为未知阻抗。采用 ph x 表示声压 试函数， 分别在内区域和外区域对式 1 进行加权 积分得到 ∫ VinnerW Δ 2ph x dV k2∫ VinnerWp h x dV 0， x ∈ Vinner 5 ∫ VouterW Δ 2ph x dV k2∫ VouterWp h x dV 0， x ∈ Vouter 6 式中，W 为加权函数。对式 5 和式 6 左边第一部分 在内、 外区域分别运用格林公式得到 在运用格林公式 时， 区域边界法向量为外法向量 ∫ VinnerW Δ 2phdV － ∫ SγW p h n γ dS － ∫ SbW p h n bdS － ∫ Vinner Δ W Δ phdV 7 ∫ VouterW Δ 2phdV ∫ S∞ W p h n ∞ dS ∫ SγW p h n bdS － ∫ Vouter Δ W Δ phdV 8 式中，n∞为无穷远处边界 S∞的外法向， 方向如图 1 所 示。将式 7 和式 8 分别代入式 5 和式 6 ， 并使用 式 2 和式 4 得到 ∫ Vinner Δ W Δ phdV － k2∫ VinnerWp hdV － ikρc∫ SbWv x dS － ikρc∫SγW ph Z dS 0 9 ∫ Vouter Δ W Δ phdV － k2∫ VouterWp hdV － ∫ S∞ W p h n ∞ dS ikρc∫ SγW ph Z dS 0 10 43振 动 与 冲 击2020 年第 39 卷 ChaoXing 2声压形函数构造 2． 1内区域声压形函数构造 在内区域， 采用径向基点插值法来构造声压形函 数， 使用添加有多项式的 RPIM 表示的声压 p x 的近 似函数为 ph x∑ n i 1 Ri x ai∑ m j 1 pj x bj RT x a pT x b 11 式中Ri x i 1， 2， ， n 为径向基函数;n 为使用 的径向基数目;pj x 为空间坐标系 xT［ x，y，z］ 中的 单项式;m 为使用的多项式数量;ai和 bj为待求常系 数;RT［ R1 x ，R2 x ， ，Rn x ］ ;pT ［ p1 x ， p2 x ， ，pm x ］ ;aT［ a1，a2， ，an］ ;bT［ b1， b2， ，bm］ 。在式 11 中， 不总是使用多项式， 当 m 0， 此时只使用基函数; 有研究表明 使用多项式能够提 高数值计算的精度和稳定性， 降低对形参数的灵敏性， 使得可以灵活广泛地选取形参数。复合二次 Multi- Quadrics，MQ ， Gaussian EXP 及薄板样条 Thin Plate Spline，TPS 函数是最常用的径向基函数 Radial Basis Function，RBF ， 如式 12 所示 MQRi x r2 i αcd 2q EXPRi x e －α c ri/d 2 TPSRi x rη { i 12 式中αc≥0， q 和 η 为形参数;d 为计算点 x 的局部支 持域中的节点平均距离;rixi－ x 为第 i 个节点 xi 与计算点 x 的距离。可见径向基函数仅是 ri的函数， 为获得良好的计算精度， 需预先确定这些形参数， 对于 一类型问题， 往往需通过数值试验来确定这些形参数 的可取值。 为求解式 11 中的 a 和 b， 需使用 n 个位于计算点 x 的局部支持域 Ωx中的场节点。支持域为一个计算点 x 常为积分点 处进行无网格插值所需要的区域， 支持 域的中心为该计算点， 为确定这组场节点， 将用到影响 域， 其被定义为一个场节点的影响区域， 该域的中心为 该场节点。当一个计算点或积分点位于某一场节点的 影响域内时， 这个场节点位于计算点的支持域内， 并被 用于构造计算点处的形函数， 图 2 显示为支持域和影 响域两者的区别。使 n 个位于计算点 x 的局部支持域 中的节点满足式 11 ， 并使用 m 个附加方程 ∑ n i 1 pj xi ai 0 j 1， 2， ， m 13 形成如下联立方程 ps [ ] 0 RPm PT m [] 0 a [ ] b G a [ ] b 14 式中pT s ［ p1， p2， ， pn］ 为声压值向量;R 及 Pm分别 为 R R1 r1 Rn r1  R1 rn Rn rn         nn 15 Pm 1 x1y1z1 pm x1  1 xnynzn pm xn         nm 16 图 2计算点 x 的支持域与节点 x i的影响域 Fig． 2Support domain of point of interest at x and influence domain of field node xi 因为对任意分布的节点， R －1通常存在并且式 11 中的多项式阶数较低， 从而在 RPIM 中不存在奇异问 题。求解式 14 中的［ a;b］ T 并将其代入式 11 得到 ph x ［ RT x ， PT x ］ G －1 ps [ ] 0 NT x ps 17 式中，NT x［ N1 x ， N2 x ， ， Nn x ］为计算点 局部支持域中的 n 个场节点的 RPIM 声压形函数。在 内区域， 式 9 中的权函数 Wi x 取为 Ni x 。 场节点的影响域大小对无网格法的计算精度影响 很大， 圆形和矩形 球形和方体型 是最为常用的二维 或三维场节点的影响域。对于圆形和球形影响域， 通 常使用影响域半径 rw来表示其大小; 对于矩形 长方 体 影响域分别在 x 和 y 方向使用 rwx和 rwy 在 x， y 和 z 方向使用 rwx， rwy和 rwz 描述其大小， 通常采用 rwx rwy 和 rwx rwy rwz来简化计算过程。r w定义为 rw α sd， 其 中 αs为影响域的无量纲尺寸; d 为节点 xi附近节点的 平均距离。由 rw定义可知， αs直接决定影响域的大小， 对于一类型问题需要通过数值试验来预先确定合适的 影响域尺寸大小。本文将采用方形或方体影响域。 2． 2外区域声压形函数构造 这一小节只对标准的 WEE 构造满足声压在无限 域自由场衰减的形函数原理进行简要概述， 关于 WEE 更为详细的知识可参阅 Cremers 等和 Moheit 等的研究。 53第 10 期吴绍维等预报振动噪声的径向基点插值无网格与无限元耦合方法 ChaoXing 结构在无限域自由场做简谐振动， 在半径为 r0恰好包 裹结构最小的球面外部区域， 声压在球坐标系中可多 极展开表示为 p r e －ikr r ∑ ∞ n 0 fn φ， θ rn 18 式中， r，φ，θ 是以 r0为原点的球坐标系。对任意的 区域 r≥r0 ε ＞ r0， 式 18 绝对且均匀的收敛， 当只取 此无穷级数的前 m 项时， 可获得 r≥r0外部区域声压近 似解。对于二维或三维问题， 标准的 WEE 是沿径向无 穷远处延伸且发散的四边形 六面体 ， 其可由位于 s- t s- v- t 坐标空间的边长为 2 正方形 正方体 母单元映 射形成， 图 3 所示为三维标准 WEE 及其映射母单元。 图 3标准 WEE 及其映射母单元 Fig． 3Standard WEE and infinite geometry mapping 沿无限边缘的逆映射为 t 1 － 2 ai r ， i 1， 2， 3， 4 19 式中r 为无限边缘上的点距离虚拟点源 1， 2， 3或 4 的距离;ai为交界面处的几何节点距相应虚拟点源的 距离。由式 19 可知， 声压 p 在母单元 t 方向轴上的 n 阶多项式形函数通过逆映射将在原单元径向方向形成 1/rn形式的展开。因此， 在母单元 t 轴上采用 n 阶多项 式可近似模拟声压幅值的衰减。n 阶的 Lagrangian 多 项式可用作构造声压幅值形函数的径向部分， 其可由 母单元 t 轴上均匀分布在几何节点之间的 n 个声节点 确定， 如图 3 所示。外区域任意一点处的声压使用所 在 WEE 形函数{ φ x } 插值确定为 ph WEE x∑ n j 1 φj x pj 20 式中φj x 为与声节点 j 对应的声压形函数， 该函数 满足 Kronecker delta 性质， 其详细表达形式可参阅文献 ［ 32］ ;pj为声节点 j 处的声压值。当计算点 x 位于 Sγ 时， { φ x } 为由几何节点 I、 II、 III、 IV 构成的有限元单 元有形函数。式 10 中的权函数 Wi x 取 G t φ s， v， t * ， 其中 G t 1 － t /2 2 为几何加权因 子， 以确保系统矩阵元素积分有限及 S∞上的积分 为零。 2． 3内区域声压形函数重构 声压由结构边界表面向外辐射的过程中是连续 的， 因此在内区域与外区域的交界面 Sγ上声压必须满 足连续和一致性要求， 即在 Sγ上满足 pinner pouter。当 计算点 x 位于 Sγ上时， 式 17 计算的 ph x 不同于式 20 计算的 ph WEE x ， 因此内区域无网格与外区域波包 络单元方程不能直接组装形成总系统矩阵方程。文献 ［ 33］ 在交界面上构造了一层交界面单元 Interface Ele- ments，IE ， 形成满足 Kronecker delta 条件的混合声压 形函数来满足声压在 Sγ上的连续和一致性， 这些交 界面单元位于沿着交界面且包含于内区域中的一层 子区域 VI中， 如图 4 所示。在所述耦合方法中， 需 在内区域为 WEE 在交界面上的几何节点配置附加 节点来形成 IE， 这使得前处理过程复杂。这些配置 的节点参与混合形函数的构造， 增加了内区域场节 点数量， 而且重构的内区域声压形函数形式复杂， 导 致计算效率下降。 图 4用于耦合 RPIM 和 WEE 的交界面单元 Fig． 4Interface elements used for coupling RPIM and WEE 为克服上述耦合缺陷， 提出一种基于改进型 WEE Improved WEE，IWEE 的混合声压形函数构造法， 该 IWEE 由标准的 WEE 和虚拟有限元网格组成， 虚拟网 格沿着交界面且位于内区域中的一层子区域 VI中与 WEE 连接， 无需在内区域为虚拟网格配置附加节点。 为述说简明， 本文采用线性虚拟网格， 如图 5 所示。虚 拟网格由位于 s- v- t 坐标空间的边长为 2 正方体母 单元映射形成。当 x 位于 Sγ时， 可使用式 17 或式 20 计算 x 处的声压， 即 ph RPIM x 或 p h WEE x 。由前述 知， 采用式 20 计算 x 处的声压时， ph WEE x 可由式 21 计算 63振 动 与 冲 击2020 年第 39 卷 ChaoXing ph FE x∑ 4 i 1 i x pi 21 图 5改进的 WEE 及其几何映射 Fig． 5Improved WEE and its geometry mapping 式中i x 为几何节点 I、 II、 III、 IV 的有限元形函数; pi为 WEE 在 Sγ上的节点声压，即 pi ph WEE xi 。当 x∈VI， 令 x0表示 x 在 Sγ上的投影。构造在 Sγ上满足 连续和一致性声压的思想是 对 x ∈ VI， 其声压由 ph RPIM x 与 p h FE x0 在虚拟网格母单元投影方向 t轴上 进行有限元插值确定， 来满足声压在 Sγ上的连续性要 求。基于这种思想， 构造内区域任意一计算点 x 处的 声压为 ph x 1 － t 2 ph RPIM x 1 t 2 ph FE x0 ，x ∈ VI ph RPIM x ， x ∈ Vinner－ VI { 22 由式 22 可知， 在计算点由内区域向外区域移动过程 中， 满足声压在交界面 Sγ上的连续和一致性要求。由 式 21 和式 22 可知， 只有 WEE 在 Sγ上的几何节点 用于构造混合声压， 这就是有限元网格被称为虚拟有 限元网格的原因。对式 22 整理得到 ph x∑ i i x pi 23 其中， i x 1 － t 2 Ni x 1 t 2 i x0 ，x ∈ VI Ni x ，x ∈ Vinner－ VI { 24 基于构造的这种声形函数， 内区域无网格与外区 域可变阶波包络单元方程可耦合组装总系统矩阵 方程。 3总系统离散方程 由于声压在交界面 Sγ上是连续一致的， 式 9 和 10 含有的未知阻尼项因大小相等符号相反而抵消， 关于内区域和外区域的系统方程祥见吴绍维等的研 究， 总系统方程可表示为 ［ K － k2M］ { p} FRPIM [] 0 25 式中K 和 M 分别为声刚度和质量矩阵;对 x∈Vinner， Kij KRPIM ij ，Mij MRPIM ij 或对 x∈Vouter， Kij KWEE ij ， Mij MWEE ij ， 其中 x 表示积分点， 下标 i 和 j 为位于 x 支持域中 的场节点 xi和 xj在全局场节点中的排序编号;FRPIM为 与结构边界节点速度相关的载荷矩阵;{ p} 为由内区 域节点声压参数或声压值组成的向量。关于 KWEE ij 和 MWEE ij 的详细计算可参阅吴绍维等的研究，FRPIM i ， KRPIM ij 及 MRPIM ij 由式 26 计算 FRPIM i ikρc∫ Sbiv x dS KRPIM ij ∫ Vinner Δ i Δ jdV，MRPIM ij ∫ VouterijdV 26 为计算 FRPIM i ， 使用一组背景网格离散结构的边界 表面 Sb。需要注意的是 背景网格实际不存在， 仅仅用 来确定积分点， 背景网格节点与问题域中的场节点独 立， 不参与场变量插值近似， 可以任意划分产生背景网 格。为采用高斯 － 勒让德求积公式进行积分， 将全局 坐标系下的不规则背景网格通过等参变换转换为局部 坐标系下的边长为 2 的标准网格 如图 6 a 和图 6 b ， 来确定高斯积分点。FRPIM i 的数值形式为 FRPIM i ikρc∑ nb l 1∫ Sliv x dS ikρc∑ nb l 1 ∑ ns h 1 ∑ nv g 1 whwgiv sh， vgdet JlJT l 槡 27 式中nb为结构边界背景网格数;ns 和 nv分别为在 s 和 v 方向上的积分点数;wh和 wg为积分点所对应的 积分因子;Jl为第 l 个背景网格转换为标准网格对应 的雅可比矩阵。 当场节点 xi和 xj的影响域之间无重叠区域时， 内 区域任何计算点或积分点的支持域都不同时包含这两 点， 从而 KRPIM ij 和 MRPIM ij 为 0。基于此规律， 采用方形或 方体影响域时， 可通过场节点影响域构造自适应背景 积分单元， 利用两场节点的影响域重叠区域 Vo作为背 景积分网格单元， 如图 7 所示， 然后对其进行等参变换 为规则的网格 如图 6 c 和图 6 d ， 采用高斯 － 勒 让德求积公式进行数值求解， KRPIM ij 的 MRPIM ij 的数值形 式为 KRPIM ij ∫ vo Δ i Δ jdV ∑ ns h 1 ∑ nv g 1 ∑ nt q 1 whwgwq 73第 10 期吴绍维等预报振动噪声的径向基点插值无网格与无限元耦合方法 ChaoXing  i s ，  i v ，  i  [] t J －1TJ－1  j s ，  j v ，  j  [] t T ‖J‖ 28 MRPIM ij ∫ voijdV ∑ ns h 1 ∑ nv g 1 ∑ nt q 1 whwgwqij‖J‖ 29 式中ns， nv和 nt分别为在 s， v 和 t 方向上的积分点 数;wh， wg和 wq为积分点所对应的积分因子;J 为背 景网格转换为标准网格对应的雅可比矩阵。 图 6整体坐标系与局部坐标系下的背景网格转换 Fig． 6Transation between global and local coordinate systems for background cells 图 7自适应背景积分单元 Fig． 7The adaptive background cell for integral 在传统的 FEM 耦合 WEE 方法中， WEE 必须与有 限元单元匹配连接， 当对有限元单元进行细化或改变 以提高声场预报精度时， 需对 WEE 进行重构并且需再 次组装系统方程， 这将增加前处理和计算时间。而在 本文提出的新方法中， 内区域场节点与外区域 WEE 无 需单元匹配连接， 由总系统方程组装可知， 外区域 WEE 的系统方程独立于内区域的系统方程; 当内区域改变 离散点分布及数量时， 无需对外区域进行单元重构， 不 必再次计算组装外区域的系统方程。 4数值算例 设计了一个圆柱壳型虚拟结构用于仿真， 其形状 和尺寸如图 8 所示， 其中系统坐标原点位于圆柱壳纵 向长度的中点处。仿真模型所处声环境参数为 声介 质密度 ρ 1． 21 kg/m3， 声速 c 343 m/s。此结构在自 由场中振动辐射的声压不存在解析解， 但是对于特定 的振速分布， 可获得准确解。在虚拟结构内部放置一 点源， 根据声压与速度的关系， 可计算出点源在虚拟结 构表面位置处的声介质质点振速， 若此虚拟结构以这 种振速分布进行振动， 其辐射的声场与点源辐射的声 场等效， 因此点源辐射的声压可作为准确值， 该点源称 为标准点源 ［34 ］。在此算例中， 将三极子标准点源置放 于坐标原点。 图 8圆柱壳结构形状与尺寸 Fig． 8Shape and sizes of the cylindrical shell RPIM 的常用 RBF 均包含有形参数， 这些形参数和 影响域的无量纲尺寸 αs共同影响该耦合方法的计算精 度。Wu 等对内区域采用 RPIM 时的收敛、 计算效率和 计算精度， αs和不同径向基函数的形参数对计算精度 的影响进行了系统研究，