STL实体的有限差分网格高效剖分算法

王跃平，薛 祥

(哈尔滨工业大学材料科学与工程学院，150001哈尔滨)

计算机数值模拟技术促进了铸造学科高速发展，为提高和确保铸件质量开创了新局面.目前，国际上在铸造过程数值模拟领域都与网格剖分前处理技术密不可分［1］.精确合理的网格剖分是提高模拟精度、减少模拟时间的重要保障.

网格剖分模块可接受的输入数据一般为各种商品化造型平台均支持的标准格式，其中STL(Stereo Lithography)格式［2-3］应用最为广泛.目前比较成熟的有限差分网格剖分算法有:切片法［4-5］、射线穿透法［6］、优化分层算法［7］等.本文在已有算法以及实现非均匀网格剖分、消除平行面误差等提高剖分精度的诸多研究［8-9］基础上，研究了铸造浇注系统有限差分网格高效精确剖分的算法，编制了程序并给出了网格剖分实例.

1 基于STL的有限差分网格剖分实现

该算法基于切片线扫描［10］原理进行有限差分非均匀网格剖分.首先，读取铸造过程数值模拟对象:铸件、铸型、浇冒口、冷铁等各STL文件作相应分析，对各实体各方向的平行面以及最小、最大坐标值排序后作为实体的剖分段边界;然后，根据实体实际情况设置并自动调整剖分步长;最后，对各STL实体的剖分段进行网格剖分.

1.1STL文件读取分析

STL文件格式将物体表面划分成很多空间小三角形，用这些三角形小平面来逼近原CAD实体，每个三角形用一个从实体内指向外部的单位矢量和3个顶点坐标表示.STL文件有ASCII与Binary两种格式，其中ASCII格式STL文件中含有“facet”关键字，而二进制文件则不具有此关键字，所以“facet”关键字可以作为判断文件格式的依据，从而采用相应的读入接口程序.本研究采用动态链表存储铸件、铸型、浇口等铸造模拟对象的STL数据为

struct STLData

{int flag［3］;//三角形面片属性标识

float facetnormal［3］;//三角形面片的法向量float vertex［3］［3］;//三角形的3个顶点坐标struct STLData*next;

}*Dhead［NoSTLs］;//各STL数据的指针数组

对STL数据进行分析，用flag值作为小三角形面片属性标识:如果3个顶点Z方向的坐标值均相等(zmax=zmid=zmin)，表明三角形为Z方向平行面，flag［2］=1;如果zmax＞zmid＞zmin，flag［2］=2;如果zmax=zmid＞zmin，flag［2］=3;如果zmax＞zmid=zmin，flag［2］=4.如果3个顶点Y方向坐标值均相等，表明该三角形为Y方向平行面，flag［1］=1;如果3个顶点X方向坐标值均相等，表明该三角形为X方向平行面，flag［0］=1.

从模拟对象的STLData动态链表结构中提取实体的X、Y、Z平行面坐标以及最大、最小坐标(Xmin、Xmax、Ymin、Ymax、Zmin、Zmax)分别按照X、Y、Z从小到大排序并剔除重复点后存储在动态链表中为

struct STLParalelPlane

{float PPlane;//平行面坐标值

int No;//某方向剖分的网格起始编号

struct STLParalelPlane*next;

}*XPPhead［NoSTLs］，*YPPhead［NoSTLs］，*ZPPhead［NoSTLs］，*AssemXPPhead，*AssemYPPhead，*AssemZPPhead;//各STL模型平行面数据(包括边界(Box［］)处退化成点或线的边界)的指针数组，以及铸造系统装配体平行面数据的指针

为了使所有的铸造模拟对象网格编号(i，j，k)统一，把铸件、铸型等各STL实体的平行面数据分别按照X、Y、Z从小到大再次排序，剔除重复点后存储在平行面链表AssemXPPhead、AssemYPPhead、AssemZPPhead中.把PPlane作为剖分段的边界与网格单元的边界，能够有效避免平行面处误差，提高剖分精度;No通过程序自动调整用户设置的网格剖分步长确定.在处理平行面数据时，可能会遇到两个PPlane非常接近的情况.经分析认为该问题是由于模拟对象的装配体文件转换成STL文件时的转换精度和算法精度不高造成的.用户可根据需要设置一误差来解决.例如设误差为0.01 mm，若两平行面数据之差小于该值，则合并为一个平行面.

程序自动调整剖分步长的基本原理为:设某剖分段的长度为Li=(PPlane)i+1-(PPlane)i，用户设置的剖分步长为dl，如果Li/dl不为整数，即表明剖分段的边界与网格单元边界不重合，需要对分区内的步长进行自动调整为

式中:f(x)为取x的四舍五入的函数;dl'为自动调整后的剖分步长.某方向上网格参数设置、调整的示例如图1所示，设用户设置的剖分步长均为dl=5 mm.

图1 某方向上网格参数示意图

网格剖分参数设置完成后，即可得各STL实体在Z、Y、X方向的起止网格编号(Nsz、Nez、Nsy、Ney、Nsx、Nex)及起始坐标(zs、ys、xs)，然后按照Z、Y、X向的顺序进行剖分.

1.2Z方向切割实体

Z向切割平面z=zk切割实体，获得一个或多个封闭轮廓线，用struct Section存储这些线段为

struct Section

{double vertex0［2］;//线段第1个端点x、y坐标

double vertex1［2］;//线段第2个端点x、y坐标

double SignofNormaly;//线段所在三角形面片法向量的Y向分量:Normaly

struct Section*next;};

切割从第n个STL(n=1，2，3，…)的最小尺寸处z=zk=zs=ZPPhead［n］+dz［k］×0.5(k=Nsz)开始，之后用z=zk=zs+(dz［k］+dz［k+1］)×0.5(k=Nsz+1，Nsz+2，…，Nez-1)平面切割实体的三角形面片.在判断三角形面片与切割平面z=zk是否相交时，先选择非Z向平行面，即struct STLData中flag［2］≠1的三角形面片，再比较三角形3顶点Z坐标值与zk值:若zmin≤zk≤zmax，则该面片与z=zk切割面相交，否则不相交.设三角形被切割棱边的顶点坐标为(x1，y1，z1)、(x2，y2，z2)，z=zk切割面与三角形棱边的交点坐标可以求得:

把三角形面片与Z向切割面交线的两个端点数据存储在vertex0、vertex1中，并把当前三角形法向量的Y向分量赋给SignofNormaly.

1.3Y轴方向切割封闭轮廓线

Y向切割线y=yj切割z=zk面上封闭轮廓线，得到与一系列成对进出的扫描线段.用struct Scanbeam记录这些线段端点的x坐标以及交点类型为

struct Scanbeam

{double vertexp;//Y扫描线与封闭轮廓线的线段交点的x坐标

int flag;//y扫描线与封闭轮廓线的线段相交点类型

Y向切割从第n个STL实体(n=1，2，3，…)的最小尺寸处y=yj=ys=YPPhead［n］+dy［j］×0.5(j=Nsy)开始，之后用y=yj=ys+(dzy［j］+dy［j+1］)×0.5(j=Nsy+1，Nsy+2，…，Ney-1)扫描线切割z=zk面上的封闭轮廓线，得到一系列交点，交点坐标可由式(1)求得.设struct Section中记录线段的两个端点的y坐标分别为ymin、ymax，且ymin≤ymax，若ymin≤yj≤ymax，则该线段与y=yj扫描线相交，将交点类型赋给flag:1-与含有ymin的端点重合;0-交点在线段中间部分;-1-与含有ymax的端点重合;2-与线段重合.

将所得交点按x坐标升序排序、处理奇异点、按奇偶配对，获得成对进出的扫描线段.

1.4X轴方向切割扫描线段

X向切割扫描线段，得到STL实体在扫描线段上的网格体心坐标(xi，yj，zk).若x=xi在成对进出的扫描线段上，则该网格为实体内部的网格.根据多材质系统网格剖分时材质属性赋值方法［11］，为材质编号为默认值(空气，通常设为0)的当前单元赋相应的材质编号值，保存网格信息，网格剖分基本完成.算法流程如图2所示.

图2 网格剖析分流程图

2 网格剖分的关键点

2.1 动态数据结构的使用

网格剖分程序运行之前，并不能够确切地知道各STL实体的三角形面片数量、需要划分的剖分段数以及XYZ各方向的网格数等信息，这就给程序设计带来了一些不便.在计算机内存资源一定的情况下，本程序采用动态链表和动态数组［12］处理网格剖分相关数据，能够较大限度地提高网格剖分数量，增强了剖分程序的稳定性和实用性.

2.2 奇异点处理

奇异点是指切割面或扫描线与STL实体相交的特殊交点，可分为:1)Z向切割STL实体时产生，如图3(a)、(c)、(d)等3种情况;2)Y向扫描线与封闭轮廓线相交时产生，如图4(b)、(c)等2种情况.

Z向切割奇异点处理为:图3(a)中轮廓退化为点的情况，直接删除;图3(c)中轮廓退化成直线或两点之间产生一条以上交线的情况，通过将切割面偏移一微小量消除;Z向平行面作为剖分段边界，可以排除情况如图3(d).通过上述处理，排除了封闭轮廓线退化成点或线的情况，避免了封闭轮廓线中出现重复线段的情况，保证了Z向切割面上的封闭轮廓线链表中记录的线段均唯一且l＞0.这样，程序更为简单，且能够有效地降低后续剖分的出错几率.

图3 切割平面与STL三角形面片相交情况

图4 扫描线与封闭轮廓线交点情况

2.3 网格剖分对STL文件的容错处理

上述剖分算法原理简单，在STL文件完好的前提下十分可靠，但可能因转换精度以及其他算法精度的原因，STL格式文件常包含各种细小的错误［13］，其中影响最大的是STL文件存在小孔洞.因此要保证程序的稳定性，要么修复STL，要么对这些错误进行容错处理.本研究采用的是容错处理方法.其基本原理为:若某条扫描线穿过封闭环上的小孔洞造成交点不配对的情况时，则给出警告:1)提示用户调整Y方向剖分步长，直到满足配对条件为止;2)保留体心坐标在该扫描线上的所有网格的原始材质(空气，材质编号为0)，待网格剖分完成后，再对网格进行过滤，纠正STL小孔洞造成的材质不匹配的网格单元.其基本原理为:针对材质编号为初始值0的网格单元，搜索并判断其3个方向:X、Y、Z上6个相邻网格单元的材质，如果某一方向上的两个相邻单元材质相同，则把该网格属性重置为相邻网格的材质编号.原理示意图如图5所示.

图5 纠正STL小孔洞导致网格材质不匹配原理示意图

3 有限差分网格剖分实例

应用上述网格剖分算法编制了三维有限差分网格剖分程序.该程序能够对各种复杂形状的STL实体进行网格剖分，剖分速度快、精度高、容错能力强.图6为剖分实例，该铸造模拟对象包括铸件、浇冒口等多个STL实体，其中，铸件轮廓尺寸为740 mm×885 mm×488 mm.图6(a)中的各实体主要由曲面组成，并且含有相对尺寸较小的倒圆、倒角、圆孔、凸台等.若采用均匀网格，并且要保留尽可能多的几何特征，就必须采用较小的剖分步长.这样网格数量会惊人地庞大，从而导致后续模拟计算时间显著增加.本算法会自动在含有小几何特征的剖分段内采用较小的剖分步长，而其他剖分段内采用较大的剖分步长，在不显著增加网格数量的前提下进行剖分，且不会出现小特征丢失或网格不连通的情况.网格剖分结果如图6(b)所示，部分参数如图6(c)所示.设输入的最大空间步长为:Δx=10 mm，Δy=10 mm，Δz=9 mm，在Pentium(R)4 CPU 1.8 GHz，内存512 MB的普通计算机上剖分耗时1.482 s，其中铸件网格数为1 007 115个.

4 结论

1)该算法能对多STL模拟对象进行分区域非均匀剖分，程序能够自动调整剖分步长、有效消除平行面处误差，提高剖分精度.

2)采用动态链表和动态数组处理数据、吸纳了现有算法处理奇异点的长处并进行改进，剖分原理简单，易于编程实现，在保证剖分精度的前提下，资源占用少，计算速度快.

3)程序在网格剖分完成后，通过网格材质匹配对STL文件的进行容错处理，能够有效修复STL小孔洞造成的扫描线进出点不配对的错误，剖分结果精确、实用性强.

［1］INOUE T，JU D Y.Analysis of solidification and viscoplastic stresses incorporating a moving boundary:an application to simulation of the casting process［J］.Journal of Thermal Stress，1992，15(1):109-128.

［2］ROSCOE L E.Sterolithography interface specification［M］.Valencia，CA:America-3D systems Inc，1988.

［3］侯华，毛红奎，张国伟.铸造过程的计算机模拟［M］.北京:国防工业出版社，2008:22-26.

［4］贾宝仟，柳百成.铸造数值模拟中直角六面体网格自动剖分技术［J］.热加工工艺，1996，(6):27-29.

［5］邱宗文，孙国雄，刘永刚，等.基于STL的切片法有限差分网格剖分的研究［J］.特种铸造技有色合金，2003，(6):28-30.

［6］周建兴，刘瑞祥，陈立亮，等.基于STL的射线穿透法网格剖分的研究［J］.铸造技术，2001，(1):15-17.

［7］徐宏，钟雪友，程军.铸件三维有限差分网格自动剖分技术［J］.铸造，2000，49(12):903-907.

［8］HUANG S H，ZHANG L C，HAN M.An effective error-tolerance slicing algorithm for STL files［J］.The International Journal of Advanced Manufacturing Technology，2002，20(5):363-367.

［9］骆国建.铸造过程计算机模拟前处理技术研究［D］.武汉:武汉理工大学，2010.

［10］梁英业，戴挺，赵建新，等.基于STL的切片线扫描法网格剖分技术［J］.铸造，2005，54(10):1002-1005.

［11］关洋.铸件充型凝固过程数值模拟前后处理技术研究［D］.沈阳:沈阳铸造研究所，2004.

［12］陈凤翔，李汪根.C++动态数组的实现与重用［J］.计算机技术与发展，2010，20(2):79-82.

［13］SZILVÁSI-NAGY M，MÁTYÁSI G Y.Analysis of STL files［J］.Mathematical and Computer Modeling，2003，38(10):945-960.