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.