基于点云数据的声学边界元分析方法

庄 超， 刘汉光， 康凯旋， 苏俊收

(1.徐工集团 江苏徐州工程机械研究院，江苏 徐州 221004；2.高端工程机械智能制造国家重点实验室，江苏 徐州 221004)

声学边界元法[1-2]，作为声学仿真的一种数值计算方法，相比于声学有限元法的技术特点在于该方法只需计算域的表面网格，这一特点对于计算大规模声学问题相对于声学有限元法具有无可比拟的优势，比如计算船体的辐射噪声，有限元法须要划分计算域的体网格，其计算量非常大，而边界元法只需获取结构表面网格即可实现辐射噪声计算，计算量大幅度下降，同时快速多级边界元法可以提高计算求解速度。针对工程广泛应用的圆柱壳结构，张俊杰等[3]探讨了边界元软件 Sysnoise两种加载方式(即基于单元(Element)加载和基于节点(Node)加载)对声辐射影响，认为基于单元加载方式得到的结果要好于基于节点加载方式，并且基于单元加载方式网格加密时收敛到解析解，而基于节点加载方式则不一定。目前，传统的声学边界元仿真过程包含建立分析对象的CAD模型，划分声学表面网格，施加边界条件并实施声学边界元数值计算等。然而，对于以下两种情况，传统声学边界元仿真呈现明显地局限性：①分析对象无现成CAD模型，须要测绘或逆向获取，无论是测绘建模还是通过逆向获取点云数据再蒙皮建立CAD模型，都将花费大量人力和时间，尤其当分析对象包含大量复杂曲面结构；②分析对象尺寸规模较大并包含复杂几何结构，声学表面网格划分将带来较大的工作量和技术难度。

逆向工程[4]，借助三维扫描和数据采集处理技术，将实物点云数据，通过三维模型重构获取三维模型，广泛应用于汽车、航空、机械工业产品开发设计，以及生物医学领域的器官三维模型重建等各领域。尤其针对具有复杂曲面的实物模型，相比于传统的测绘方法，具有天然的效率和优势。基于点云数据，技术人员将逆向工程与数值分析相结合，前者获取的三维网格模型作为后者的计算网格数据，实现两者的无缝衔接。例如，在生物医学应用领域，赵岩等[5-6]，将CT扫描图像在Mimics等逆向工程软件进行三维模型重建，并应用Remesh优化三维面网格以获取规则化的三角面片，然后导入有限元软件中划分体网格，赋以材料属性并添加边界条件，实现有限元数值分析。

考虑到传统声学边界元的上述缺陷，以及逆向工程与数值分析相结合的可行性，本文首次将逆向工程和声学边界元法相结合，将前者获取的点云数据转换为可用于计算的声学边界元表面网格，实现两者的无缝衔接；同时，给出了小型扫路机厢体声传播问题的具体算例，并与传统声学边界元仿真结果对比，得到了较高的吻合度，证明了本文方法的有效性与可行性。

1 点云及三角网格获取

利用三维激光扫描仪获取小型扫路机厢体部件(见图1)的点云数据，在Geomagic逆向处理软件中，同实物比对以获取连续的、无离散数据的点云为原则；通过直接观察法、弦高差法，采用自动修复、去除特征、清除、填充孔等操作，识别并去除点云的噪声点和冗余数据[7-8]。

图1 小型扫路机Fig.1 Small sweeper

考虑计算量基于均匀采样或者弦偏差采样实现点云精简，采用曲线、曲面插值补充法对缺失的几何特征完成数据补齐；基于曲率均匀法借助快速光顺对其进行滤波光顺处理[9]，获得最终的点云数据三角网格，为后续的数值计算提供最基本的网格数据，该三角网格是连续的且能够表达实物几何特征，如图2所示。图2中点云三角化模型为了保留实体几何特征，曲面网格数量较大，整个模型范围内网格单元分布不均与、网格密度相差大，同时三角面片大小与形状不均匀，存在大量狭长型、内角多为小锐角和大钝角的三角面片，以及阶数大于6的顶点。

图2 三角网格模型Fig.2 Triangular mesh model

2 三角网格的重网格化

由图2可知，三角网格的质量并不满足插值与积分计算的需求，无法直接用于边界元数值计算，须要进行网格精简和单元质量优化处理，以获取既符合几何表达又满足数值计算要求的三角网格单元模型。沈建国等[10]提出面向有限元分析的三角网格迭代优化，通过网格细分、简化和规则化处理，在给定的误差范围内有效地改善三角形形态，提高了网格质量，为有限元单元的生成提供合适的网格模型。

为实现三角网格的重网格化[11]；首先，根据数值计算的类型、精度、模型整体规模等，给定目标数量或单元尺寸简化模型；然后，通过重网格化实现单元质量优化，包括几何和拓扑优化两个方面，前者以三角面片形状和尺寸均匀为目标，后者以网格顶点具有相同的阶数为目的。

三角网格的重网格化的基本思想是对于给定的初始网格，基于采样均匀性、规则化、尺寸和单元形状等准则，重新生成网格在较好近似原网格的基础上获取较高的网格单元质量。本文采用半边数据结构(Halfedge Data)以及相应的库函数组成的表示和管理多边形网格的数据结构OpenMesh[12-13]，来存储、管理以及遍历模型的网格数据，并根据Botsch等[14]提出的各向同性重网格化实现单元质量优化。其算法实现过程为：①分裂所有长度大于目标边长的边；②消除所有长度短于目标边长的边；③翻转边以获取最优的平均顶点阶数；④在切平面中移动顶点至重心位置。

将上述步骤迭代数次，获得目标边长、顶点阶数为6的三角面片网格。然后，采用基于面积的切向光滑处理，将顶点移动至重心加权的形心位置，加权系数计算如式(1)，最终得到理想的三角网格模型，具有网格单元密度与尺寸均匀、形状均接近等边三角形、绝大多数顶点阶数为6等特点。原厢体三角面片网格的重网格化后的网格模型，作为声学边界元分析的单元模型，如图3所示。

(1)

图3 三角网格的重网格化Fig.3 Remeshing of triangular mesh

3 编写数值计算文件

经重网格化获取的厢体三角网格单元模型，仍然是一种表达几何的模型，不能直接用于计算。表达几何的三维网格数据包含网格顶点的三维坐标、编号以及三角面片的拓扑连接关系等。

本文通过编程，将上述数据写成通用的、包含数值计算信息的INP格式数据文件，包含节点、单元、单元类型等基本数据。图4给出了表达几何及数值计算的数据文件格式。通过编写数值计算数据文件，我们获取了可以导入声学边界元计算程序的单元模型数据。

图4 几何与计算数据文件格式Fig.4 File formats for geometric and computational data

4 声学边界元理论

声学问题控制方程Helmholtz微分方程

2p+k2p=0

(2)

式中:k为波数，k=ω/c。控制方程包含Dirichlet,Neumann以及Robin等三类边界条件。对于稳态声场外问题，还须要满足远场Sommerfeld辐射边界条件。

(3)

借助Helmholtz方程基本解u*与格林公式[15]，

(4)

将式(2)和式(4)代入式(5)，得到间接边界元网格上的单层势σ与双层势μ表示的Helmholtz间接边界积分方程

(6)

式中:单层势σ代表间接边界元两侧的速度差即

(7)

双层势μ代表间接边界元两个的声压差即

μ=p+-p-

(8)

将式(6)中的源点P移到求解域的边界，获得边界上的声压及其外法向质点振速之间的边界积分方程。

(9)

(10)

为了数值求解边界积分方程，须要将声学域的边界划分成有限个单元，这里用二维三角形面单元。将边界离散成Ne个单元，则边界积分式(9)变成

(11)

式中：Γe为单元e的边界。对于边界上所有N个节点，可以得到N个方程，用矩阵形式表示为

(12)

式中：B,C和D为系数矩阵;σ为单层势列向量与μ为双层势列向量;f和g为外界激励向量。求解上述方程得到边界单元各节点的单层势σ与双层势μ，然后声场区域内任一点声压为

p={Aσ}T{σ}+{Aμ}T{μ}

本文用于离散计算域的二维三角形面单元，即是点云数据经重网格化获取的分布均匀、形状规则的二维三角形网格单元。

5 算 例

为了验证本文所述方法的可行性，将其用于模拟小扫路机厢体内放置点声源的声传播问题，并与传统声学边界元方法的结果进行对比。计算流体域为空气，边界条件为位于厢体内某一位置的点声源，其频谱特性如图5所示。频率范围16～1 800 Hz，在厢体左、右侧7.5 m某位置处放置测点，拾取声压响应。

首先，建立基于点云数据的声学边界元计算模型，将上述重网格化后的小扫路机厢体三角形网格模型经数据格式转换得到的INP数据模型，并用于声学边界元数值计算。

图5 点声源频谱曲线Fig.5 Spectrum curve of point sound source

然后，建立传统声学边界元分析模型。①将逆向获取的点云数据导入CAD软件Pro Enginner，通过表面蒙皮、光顺等实现三维几何重建以获得厢体的三维CAD模型，如图6所示。②利用网格划分工具HyperWorks得到用以计算的表面网格，导出NASTRAN类型的Standard format格式文件，为了满足计算精度需求，须要花费一定的时间，基于质量评价准则，通过单元质量优化来获取整体单元质量较好的网格单元模型，如图7所示。③导入声学边界元程序LMS Virtual. Lab Acoustics，定义为声学网格，定义声学材料为空气，将声学材料定义到声学网格上，定义场点，定义点声源并导入图5所示的频谱曲线，先后计算声场分布、场点响应以及场点上的声压频率响应函数。

用于仿真的计算机处理器为Intel 3.60 GHz，内存16 GB。表1给出两种声学仿真方法各步骤及时间花费，实物逆向时间相同，同时由于目标计算网格数量差不多，所以分析计算时间基本一致；主要时间差异来源于从点云数据获取计算网格的时间，传统声学边界元须建立CAD模型并划分网格，共须240 min，而基于点云数据的声学边界元方法通过程序实现点云三角网格到计算网格的转换，仅须15 min。

表1 基于点云数据的与传统声学BEM仿真时间Tab.1 Simulation time of based on point cloud data and traditional acoustic BEM min

图6 CAD模型Fig.6 CAD model

图7 边界元表面网格Fig.7 Surface mesh for boundary element

图8和图9，分别给出了上述应用基于点云数据的声学边界元与传统声学边界元仿真计算结果，即厢体左、右侧测点的声压响应频谱曲线。由图可知，两种方法仿真结果吻合度较高，这也证明了本文提出的基于点云数据的声学边界元分析方法的有效性。

图8 左侧点声压响应频谱Fig.8 Sound pressure response spectrum of the left side

图9 右侧点声压响应频谱Fig.9 Sound pressure response spectrum of the right side

6 结 论

本文提出的基于点云数据的声学边界元分析方法，将逆向工程与声学边界元数值计算相结合，把扫描获取的不规则的几何网格模型优化处理，得到形状规则的网格模型以满足数值计算的需求，然后转换为包含单元类型和拓扑关系的计算数据模型以用于声学边界元分析，实现了CAD几何模型与CAE仿真分析的无缝衔接。

于是，对于实物模型的声学分析，在通过逆向扫描获取点云数据之后，不再需要根据点云建立CAD几何模型并划分声学边界元表面网格的传统路径，只须将点云数据进行优化处理并添加单元数值计算信息，即可得到适用于仿真计算的网格单元数据，该方法对于具有复杂曲面的实物模型尤其能体现出巨大优势。文中给出的小型扫路机厢体声传播问题具体算例，与传统声学边界元仿真对比，得到了较高的吻合度，证明了该方法的可行性与有效性。