面向无网格法的质点生成算法

周词扬 郑澎 马智博 刘敏娟

摘要：

介绍一种面向无网格数值模拟方法的质点生成算法。将四面体、三角形网格生成算法分别用于空间平面、曲面和实体模型，将网格单元的属性赋予单元内某一点作为质点，并生成对应质点集。为研究不同网格生成算法和质点生成算法对质点集的影响，提出一种质量评价标准，开展对不同算法组合的质量分析，得到网格生成算法和质点生成算法中的最佳组合。

关键词：

无网格法； 质点生成； 质点集； 网格生成； 评价标准

中图分类号： TP391.9

文献标志码： B

Particle generation algorithm for meshless method

ZHOU Ciyang1a， ZHENG Peng1， MA Zhibo2， LIU Minjuan1a

（1. a. Computer Application Institute， Mianyang 621900， Sichuan， China； b. Software Center for High Performance

Numerical Simulation， Beijing 100088， China， China Academy of Engineering Physics； 2. Beijing Institute of Applied

Physics and Computational Mathematics， Beijing 100094， China）

Abstract：

A particle generation method is introduced for meshless numerical simulation method. Etrahedron and triangle mesh generation algorithm is used for plane， surface and entity model； the properties of the mesh unit are given to a point inside the unit as particle and the corresponding particle sets are generated. In order to study the influence of different mesh generation and particle generation algorithms on particle sets， a quality evaluation standard is proposed. The quality analysis of different algorithm combinations is carried out， and the optimal combination of mesh generation algorithm and particle generation algorithm is obtained.

Key words：

meshless method； particle generation； particles set； mesh generation； evaluation standard

收稿日期： 2017-12-19

修回日期： 2018-01-23

基金項目： 国防基础科研计划（C1520110002）；中国工程物理研究院双百基金（HT1609-01-KY124-019）

作者简介：

周词扬（1993—），男，湖北仙桃人，硕士研究生，研究方向为高性能数值模拟前处理，（E-mail）1095798605@qq.com

0 引 言

随着高性能计算技术的发展，有限元法在核物理学、爆轰、材料学等领域发挥着越来越大的作用。但是，有限元法在处理材料变形、高速碰撞等问题时，由于网格发生形变，重新划分困难，导致计算速度下降，计算精度降低，严重时甚至出现计算崩溃的情况。采用无网格法解决这类问题时，由于不依赖网格或只使用部分网格，可避免网格重新划分等问题，能有效解决有限元法所遇到的困难。然而，目前无网格法仍处于研究初级阶段，其相应的工程软件较少，部分功能尚未实现。[1]无网格法的基本思想是使用一系列无网格节点排列，采用一种与权函数相关的近似，使某个域上的节点可以影响研究对象任意一点的力学特性，摆脱不连续性的束缚，保证求解的精度。[2]因此，在前处理过程中如何布置和划分用于表达连续体信息的质点成为重要研究方向。本文设计一个面向无网格法、用于实现几何模型质点集生成的算法，生成模型的四面体/三角形网格，由网格得到对应的质点集。讨论质点生成算法的不同实现方法并比较其对结果的影响，通过对同一个模型进行测试，比较不同实现方案结果的差异。

采用无网格法中的光滑粒子流体动力学（smoothed partical hydrodynamics， SPH）方法，通过四面体网格得到模型质点。SPH方法的主要思想是任何一个连续系统可离散为一系列任意分布的质点，所有有关这一系统的量都被认为集中在这些质点上。[3-4]生成四面体网格后，用网格单元内的某点代替网格，网格单元的属性值集中于该点，使生成的质点集能用于SPH方法及其衍生方法的计算。

1 三维实体内部质点集生成算法

1.1 算法概述

用无网格数值模拟方法对模型进行仿真计算时，需要将模型离散为若干质点，用质点代表周围区域，质点的属性即为周围区域的属性。

生成模型的四面体网格，每个网格单元对应一个质点，质点质量为网格单元质量。此外，算法需要一个标准量化网格质点分布状况，以便于比较不同生成算法之间的效果。

质点生成算法流程为：

（1）生成模型的初始网格。

（2）细化初始网格到用户指定程度。

（3）采用质点生成算法求出细化网格每个网格单元的质心，得到模型对应的质点集。

（4）对质点集采用相应评价标准得到具体数值，与设置的阈值比较，如果不符合要求则重复步骤（2）和（3），直到满足要求为止。

在上述步骤中，网格生成算法和质点生成算法均会对最终结果产生影响，不同的算法组合会得到不同的结果。本文讨论每个算法的多种实现方式，比较不同方法所生成质点集的质量并确定一种最佳方法。

1.2 网格生成算法

本算法对四面体网格生成算法的要求是网格质量良好、网格密度可控：网格质量良好则四面体单元更接近于正四面体，相邻质点之间的距离方差更小，质点分布更加均匀；网格密度可控是为了保证能细化网格到指定程度。

网格生成算法采用Delaunay算法保证网格质量。四面体Delaunay网格的特性是对任意网格单元，其外接球内不存在除了构成网格单元顶点之外的点[5]，网格整体质量良好，畸形网格数量较少，并且能保证限定条件在网格中的一致性，网格尺度和质量可控。采用自适应网格生成算法，通过限制网格单元的体积或相邻点距离实现功能。

限制期望距离的Delaunay细化算法需要先指定一个距离L，使生成的网格边长在该值附近，然后为算法设置一个点生成规则，检查每个网格单元，如果满足点生成规则时则生成新的点并细化网格，直到所有单元均无法再细化为止。对应的点生成规则[6]为：

（1）当某条限定边不存在于网格中或网格中某条边被干涉，则将该边的中点加入到四面体网格中并重构网格。

（2）如果某个限定面不存在于网格中或网格中的三角面片被干涉，则将该面的外心加入到四面体网格中并重构网格。如果出现规则（1）中的情况，则用R1的方法处理干涉边。

（3）如果四面体的Radius-edge比大于阈值，或者外接球半径大于L，则将四面体单元的外心加入到网格中并重构网格，如果网格出现规则（1）和（2）中的情况，则分别按第1.1节中的（1）和（2）处理对应线面。

在点生成规则中，三角面片被干涉是指以该三角面片外心为球心，外接圆半径为球半径的球含有除三角形顶点外的其他点。线段被干涉是指以线段中点为球心，线段长度的1/2为半径的球含有除线段端点以外的点。

基于上述点生成规则，网格细化算法流程如下：

（1）将初始四面体网格的网格单元存储在队列中；

（2）从队列取出一个四面体单元，按照点生成规则生成新的点并重构网格，删除队列中已不存在于网格中的四面体单元，将新生成的四面体单元加入队列；（3）重复步骤（2），直到队列为空。

限制最大体积的自适应Delaunay算法需要指定一个体积V，生成的网格中所有网格单元体积均小于V。对于限制最大体积的网格生成算法，对应的點生成规则与限制距离的规则大致相同，唯一的不同之处是点生成规则中第（3）条加入外心的条件是四面体的Radius-edge比大于阈值或者体积大于V。算法流程与限制距离的完全一样。

1.3 质点生成算法

质点生成算法要求确保生成的质点在网格内。若某个网格单元对应的质点有在单元外的可能性，则对于位于边界的单元会出现对应质点在几何模型外的情况，同时内部质点分布不均匀。设计2种质点生成方法，取四面体内心和取四面体重心。

取四面体内心法：已知四面体4顶点坐标，求内心坐标。设4顶点坐标分别为A1（x1，y1，z1），A2（x2，y2，z2），A3（x3，y3，z3），A4（x4，y4，z4），内心坐标设为（xa，ya，za），则求解公式[7]为

xa=4i=1（sixi）/4i=1xi

ya=4i=1（siyi）/4i=1yi

za=4i=1（sizi）/4i=1zi

（1）

式中：si为点Ai所对侧面的面积。

取四面体重心法：已知四面体4顶点坐标，求重心坐标。设重心坐标为（xg，yg，zg），则求解公式为

xg=4i=1xi/4

yg=4i=1yi/4

zg=4i=1zi/4

（2）

1.4 质点集质量评价标准

对于生成的质点集，质点分布越均匀，后续计算中计算精度和计算效率越高，计算出错的可能性越低。同时，无网格法用质点代表对应的四面体，每个质点的质量为对应四面体的质量，质点的质量值越集中，说明模型划分得越均匀，有助于提高后续计算精度。假设模型密度恒定，则可以用体积值代替质量。因此，评价质点集的质量可从质点分布质量和质点所在单元体积两方面进行。

为了量化质点集分布质量，生成质点集的Delaunay网格，将Delaunay网格质量作为质点集质量。由于Delaunay网格与Voronoi图为对偶关系，Voronoi图中每个空间单元均为对应点的邻域，而邻域共面的两点在Delaunay网格中是相连的，因此Delaunay网格中线段两端点为空间相邻点，因此该网格的质量可表示质点集的分布质量。

网格的质量可以用能量函数[8]、平均曲率[9]等方式表示，本文用网格单元的平均质量表示网格整体质量。每个网格单元的质量评价标准如下。

Radius-edge ratio标准公式[10]为

d=r/l （3）

式中：r为四面体外接球半径；l为四面体最短边长。

ParaQ标准公式[11]为

Q=CdV/（1i

式中：V为四面体体积；lij为端点是顶点i、j的线段长度；Cd为使正四面体的质量度量值取1而采用的比例因数，取Cd=1 832.820 8。

ParaV标准公式[11]为

Q=723V/（1i

l2ij）1.5

（5）

式中：V为四面体体积；lij为端点是顶点i、j的线段长度。

每个四面体单元对应质点的属性为该四面体单元体积，单元体积越集中，属性值越集中，后续计算的精度和效率越高。用最大体积与最小体积之比衡量质点集属性的分布状况。体积比计算公式为

σ=Vmax/Vmin （6）

式中：σ为体积比；Vmax为最大四面体单元体积；Vmin为最小四面体单元体积。

得到质点集后，综合对比3个网格评价标准，比较不同方法生成的质点集的分布质量，并讨论质点集的体积比。

2 空间面片质点生成

在无网格数值模拟方法中，常常需要离散空间平面或曲面。得到的质点集一般用于作为边界条件求解碰撞发生时间、状态等；除此之外，还可能用于求解该平面或曲面的物理属性，例如求解材料变形时表面材料的状态和形状。对于前者，离散后的质点必须在平面或曲面上，质点没有属性；对于后者，质点与平面或曲面之间允许有一定距离，质点可看作将所在微元面集中于该点，该点属性代表所在微元面属性。

空间面片质点生成算法输入为一组空间平面或曲面，输出为一组质点集。空间平面的质点生成算法的流程为：

（1）生成空间平面初始网格，将这些初始网格单元按所属平面分类，以队列的形式存储。

（2）遍历初始网格队列，细化每一个外表面网格，采用质点生成方法生成每个网格单元质点并合并为质点集。

（3）重复步骤（2）直到遍历完队列所有单元，整合所有质点集为一个总质点集并输出。

空间曲面的质点生成算法流程为：

（1）生成曲面的初始三角形网格。

（2）采用曲面逼近算法生成新点并细化网格

（3）根据实际应用的不同取三角形网格的顶点、所有网格单元内心或重心的集合作为质点集。

曲面质点生成算法根据选择顶点方式不同又分为映射法和直接法2类，本文采用直接法中的网格前沿法。

2.1 平面三角形网格生成算法

由于在二维网格中Delaunay三角形网格具有空外接圆性和最优性，因此采用Delaunay算法作为三角形网格生成算法。细化网格采用限制网格边长的方法确保能细化到指定程度。

细化算法质点生成规则为：

（1）如果边被干涉或限定边不存在于网格中，则将边的中点加入网格并重新划分。

（2）如果三角形被干涉或三角形外接球半径大于指定值L，则将三角形的外心加入网格并重新划分；如果外心与某一条边干涉，则用规则（1）的方式处理干涉边。

细化算法检查每个三角形单元，如果满足上述规则，那么细化该单元。整个生成算法流程与四面体的相同，只是处理的网格單元为三角形单元。

2.2 曲面三角形网格生成算法

根据网格生成方法的不同，空间曲面网格生成算法可分为直接法和映射法。其中直接法又可以分为曲面分解法和网格前沿法2类。本文采用网格前沿法。网格前沿法[12]最早由LHNER提出，用于生成二维或三维非结构网格，其中生成三角形网格时要求面为平面。

将网格前沿法用于空间曲面三角形网格生成时，由于网格与模型不重合，因此需要对算法进行修改。LAU等[13]提出一种基于网格前沿法的空间曲面网格生成算法，该算法与平面网格不同之处是生成新边界点的方式不同。

假设边AB为下一个从前沿中删除的边，h为新网格单元预期高度，d为前沿中某条边与曲面相切且与边夹角为90°的向量，不同前沿点的d值不同，则新点坐标为

Cz=M+hd （7）

式中：M为前沿边队列中某一个端点。选择的端点不同，得到的新点不同。比较不同点与边AB构成的三角形的质量，选择质量最好的网格单元作为新网格单元，同时更新网格前沿和网格单元的队列。

2.3 质点生成算法

质点生成算法输入一个三角形网格，输出一组质点。一般情况下，每个网格单元对应一个质点，计算时该质点代表对应单元。质点生成方法与三维质点生成算法相似：如果质点有属性，则取三角形的内心或三角形的重心；如果质点没有属性值，则可以取三角形网格的顶点。

已知三角形三顶点坐标为A1（x1，y1，z1）、A2（x2，y2，z2）、A3（x3，y3，z3），设内心的坐标为Aa（xa，ya，za），重心的坐标为Ag（xg，yg，zg），则坐标值为

xg=3i=1xi/3

yg=3i=1yi/3

zg=3i=1zi/3

xa=（ax1+bx2+cx3）/（a+b+c）

ya=（ay1+by2+cy3）/（a+b+c）

za=（az1+bz2+cz3）/（a+b+c）

（8）

式中：a的值为|BC|；b的值为|AC|；c的值为|AB|。

2.4 质点集质量评价标准

二维质点集的评价标准与三维类似，不同之处是用于质量评价的网格是由模型表面生成的三角形网格，网格质量为所有网格单元的平均质量，将该值作为质点集的质量。该方法得到的质量能用于比较不同网格生成方法得到质点集之间分布质量，无法用于比较不同质点生成方式得到质点集的质量。网格单元的质量有以下几种标准。

Radius-edge ratio标准公式为

d=R/l （9）

式中：R为三角形外接圆半径；l为三角形最短边长。

纵横比标准公式为

d=R/r （10）

式中：r为三角形内切圆半径。

比较不同质点集质量时，使用上述标准得到质点集的多个质量值，综合比较这些数值，判断质点集质量的优劣。

3 实验验证

3.1 三维实体质点生成算法

爆炸装置模型见图1，该模型共5个部件，每个部件的内容物不同。无网格算法模拟装置的爆炸过程，判断装置从最上方起爆，爆轰波能否到达最下面。设半圆柱最上方部件的顶面和侧面为固壁。

采用限制最大距离的自适应Delaunay算法时，期望距离设为0.000 23 m，模型质点生成结果见图2。采用限制最大体积的自适应Delaunay算法时，设置最大体积为1 mm3，得到的质点集见图3。

图 1 爆炸装置模型

图 2 限制最大距离的自适应算法模型质点集

图 3

限制最大体积的自适应算法模型质点集

2×2×12长方体质点集生成质量见表1，按照上文介绍的3种质量评价标准进行质量分析。由此可知：内心生成的质点集分布质量优于重心的；质点集生成的网格中，四面体单元整体质量优于后者，但部分四面体单元质量较差，从而导致在用Radius-edge标准判断质量时数值过大。

表 1 2×2×12长方体质点集生成质量

整理半圆柱模型最上方部件网格的各网格单元体积分布，见图4。由此可知，原四面体各网格单元体积取值范围很大，最大值为1.750E-11，最小值为1.849E-13，最大值约为最小值的95倍。说明四面体网格生成算法得到的网格单元的大小不均匀，需要改进。

图 4 四面体网格单元体积分布

3.2 表面质点生成

生成固壁对应的质点集时，质点取对应三角形单元的重心，固壁质点见图5。

图 5 固壁质点

生成整个模型的表面质点，质点取网格单元的节点，见图6。

图 6 整个模型的表面质点

对固壁的2个表面进行质量分析，结果见表2。由此可知：由于纵横比接近于2，因此固壁上的三角形网格中每个网格单元形状接近正三角形，网格质量良好。

表 2 模型表面质点生成质量

固壁对应的三角形网格单元面积见图7。由此可知，三角形面积的最大值为7.22E-8，最小值为3.00E-8，最大值约为最小值的2.4倍。结果表明：对固壁进行三角形划分时得到的网格单元面积值集中，有助于后续计算保持高精度。

图 7 固壁对应的三角形网格单元面积

4 结束语

对多个模型生成对应的内部质点集和表面质点集。由分析结果可知：网格生成算法为限制距离的自适应Delaunay算法时，生成的质点分布质量优于限制最大体积的自适应Delaunay算法；质点生成方式中，取内心的方式得到的质点集整体分布比取重心方式的更加均匀，但局部地区分布状况不如后者。目前的算法得到的质点集仍存在问题，当给每个质点赋予对应四面体体积值时，数值分布范围过大，数值分布比较发散，采用适当的优化算法优化由模型得到的四面体网格，使网格单元体积集中于某个值附近。表面质点集的质点分布均匀，属性值集中，能用于实际SPH无网格法的计算中。

参考文献：

[1] 王宇新，孙明，张建臣.无网格MPM法三维前处理系统设计[J]. 计算力学学报， 2008， 25（3）： 392-396.

[2] 刘天祥，刘更，朱均，等.无网格法的研究进展[J]. 机械工程学报， 2002， 38（5）： 7-12.

[3] 王宇新，陈震，张洪武，等.多层抗暴结构冲击响应无网格MPM法分析[J]. 工程力学， 2007， 24（12）： 186-192.

[4] 顾元通，丁桦.无网格法及其最新进展[J]. 力学

进展， 2005， 35（3）： 323-337.

[5] 武晓波，王世新，肖春生.Delaunay三角网的生成算法研究[J]. 测绘学报， 1999， 28（1）： 28-35.

[6] SI H. On refinement of constrained Delaunay Tetrahedralizations[C]//Proceedings of the 15th International Meshing Roundtable Conference. Birmingham， USA， 2006： 509-528.

[7] 王曉凤，李永利.四面体的内心坐标公式及其应用[J].平顶山工学院学报， 2004， 13（4）： 75-78.

[8] HOPPE H， DEROSE T， DUCHAMP T， et al. Mesh optimization[J]. Conference on Computer Graphics & Interactive Techniques， 1993， 27： 19-26.

[9] KARBACHER S， HAEUSLER G. New approach for the modeling and smoothing of scattered 3D data[R]. San Jose， USA： International Society for Optical Engineering， 1998.

[10] 李海生.Delaunay三角剖分理论及可视化应用研究[M]. 哈尔滨： 哈尔滨工业大学出版社， 2010： 66-67.

[11] 聂春戈，刘剑飞，孙树立.四面体网格质量度量准则的研究[J]. 计算力学学报， 2003， 20（5）： 579-582.

[12] LHNER R， PARIKH P. Generation of three-dimensional unstructured grids by the advancing-front method[J]. International Journal for Numerical Methods in Fluids， 1988， 8（10）： 1135-1149.

[13] LAU T S， LO S H. Finite element mesh generation over analytical curved surfaces[J].Computer & Structures， 1996， 59（2）： 301-309.