劉尚蔚,胡雨晨,魏魯婷,魏 群

(華北水利水電大學 水利學院，河南 鄭州 450045)

水工建筑物中存在很多形狀特殊、結構復雜的部分，3D打印技術可以將其快速打印成型，并應用在水利施工與模型試驗等方面。隨著BIM技術的快速發(fā)展和在水利行業(yè)的廣泛應用，BIM三維模型與快速成型技術也已緊密結合。STL模型具有良好的幾何適應性，以“Polygan Soup”的方法可精確表示各種復雜結構[1]，對于實體數(shù)字圖形產(chǎn)生閉合的STL網(wǎng)片(稱為實體類STL網(wǎng)格)，通過切片技術獲得3D打印的數(shù)據(jù)格式，目前已廣泛應用于快速成型技術。這種STL三角面網(wǎng)格由于其閉合特性，包括了內(nèi)部的孔洞或凹凸體造型，有邊界、有厚度，并有不同材料屬性制約的STL網(wǎng)格集合體的網(wǎng)格錯綜復雜，無論用于是三維打印數(shù)據(jù)還是轉化為四面體或六面體單元，對模型精細化要求高，必須切割成規(guī)則的形狀，費時費力[2]。這使得利用三維數(shù)字模型進行3D打印成型并進行有限元網(wǎng)格轉化和計算，常常需要借助專業(yè)軟件進行處理。

目前常用的攝影成像或激光掃描所得到的三維模型僅僅是所見表面的模型[3]，其形成的STL網(wǎng)格是單層表面類STL網(wǎng)格，不具有厚度尺寸，這與三維建模形成的STL實體類閉合網(wǎng)片有著很大的區(qū)別。魏魯雙[4-5]等對于單層表面STL網(wǎng)格進行等厚度處理，形成等厚度3D打印模型，并轉化為有限元網(wǎng)格，將3D打印技術與有限元網(wǎng)格融合，提出了薄壁結構和鋼桁架的STL文件與轉換方法，提供了3D打印與數(shù)值計算網(wǎng)格融合一體的快速實用技術。在此基礎上，本文對非均勻厚度的表面類STL網(wǎng)格的打印數(shù)據(jù)處理提供了新方法，并以某水電站變厚度尾水管為例，得到實際驗證。

1 STL文件與有限元網(wǎng)格轉化

1.1 STL文件

STL文件格式簡單，只包含三角面片網(wǎng)格的三個角點和法矢量信息，無拓撲關系，文件所需存儲空間小，成型快，已經(jīng)成為3D打印的標準打印格式，大多數(shù)三維建模軟件都可以輸出，如AutoCAD、3DMAX、Catia、Pro/E等[6]。

圖1 STL三角面片外法線與角點圖

solid Object 01

facet normal (x,y,z)/*第一個三角單元的法矢量*/

outer loop

vertex(xi,yi,zi)/*第一個單元的第一點坐標*/

vertex(xj,yj,zj)/*第一個單元的第二點坐標*/

vertex(xk,yk,zk)/*第一個單元的第三點坐標*/

endloop

endfacet

......

facet normal(x,y,z)/*第m個三角單元的法矢量*/

outer loop

vertex(xi,yi,zi)/*第m個單元的第一點坐標*/

vertex(xj,yj,zj)/*第m個單元的第二點坐標*/

vertex(xk,yk,zk)/*第m個單元的第三點坐標*/

endloop

endfacet

endsolid Object 01

其中，“solid”表示一個模型，“Object01”則表示該模型的名稱，“endsolid”表示整個文件的結束，“facet normal”到“endfacet”間包含一個三角形面片的單元信息，“facet normal”表示三角形外法線方向單位矢量，“outer loop”到“endloop”間包含的是三角形面片三個角點的三維坐標，“endfacet”表示一個三角面片信息的結束。

1.2 STL文件的處理

STL文件在存儲數(shù)據(jù)時依次存儲每個單元的法矢量和三個角點的坐標信息，這會使得同一節(jié)點被重復存儲[10]。在轉化有限元數(shù)據(jù)格式的過程中，設置閾值(r)，將閾值范圍內(nèi)的點合并歸類，清除奇異點，處理后的點重新排列生成整體節(jié)點坐標。

某STL文件中包含m個三角面片單元，用整體單元矩陣表示如下：

(1)

根據(jù)數(shù)據(jù)精度給定閾值(r)，對閾值內(nèi)點進行合并歸類消除奇異點后，節(jié)點數(shù)由3*m個變?yōu)閚個，n個整體節(jié)點矩陣如下：

(2)

由m個單元表示整體單元矩陣如下(a、b、c、...、n1、n2、n3均為上述n個整體節(jié)點矩陣中的元素)：

(3)

1.3 STL三角網(wǎng)格轉化為有限元網(wǎng)格

STL三角網(wǎng)格作為有限元計算網(wǎng)格仍然不理想，需通過編寫程序，處理鄰近點、狹長單元，得到符合轉化要求的STL外表面網(wǎng)格，再拉伸生成有限元體網(wǎng)格。STL外表面網(wǎng)格生成及轉化為有限元網(wǎng)格的流程見圖2。

圖2 STL網(wǎng)格轉化有限元網(wǎng)格流程圖

1.3.1 三角網(wǎng)格優(yōu)化處理

STL文件中的三角網(wǎng)格存在節(jié)點重疊、間距較近、網(wǎng)格角度過大或過小等問題，造成外表面三角網(wǎng)格無法直接用于有限元體網(wǎng)格的生成，因此需要通過程序?qū)W(wǎng)格進行清除奇異點處理，以滿足有限元計算要求。

如果相鄰兩單元節(jié)點重疊，則會造成節(jié)點和單元數(shù)據(jù)重復；如果間距較近，則出現(xiàn)狹長單元等問題。給定一閾值(r)判定節(jié)點是否重疊或間距較近，節(jié)點P3、P4間距若等于0(見圖3)，則判定為重疊節(jié)點，重復節(jié)點和單元應刪除；若小于r(見圖4)，則判定為節(jié)點間距較近，需合并節(jié)點并刪除單元。合并、刪除節(jié)點和單元后將修改節(jié)點與單元編號，得到符合有限元計算，形成質(zhì)量良好的外表面網(wǎng)格(見圖5)。

圖3 兩節(jié)點重疊的三角網(wǎng)格圖

圖4 兩節(jié)點距離較近的三角網(wǎng)格圖

圖5 優(yōu)化后的三角網(wǎng)格圖

1.3.2 生成有限元網(wǎng)格

有限元外表面網(wǎng)格的節(jié)點沿著某一節(jié)點所連接的所有三角單元的合矢量的負方向拉伸，新節(jié)點構成了體單元內(nèi)表面網(wǎng)格，將歸入拉伸前的總節(jié)點編號及相應單元中。由STL網(wǎng)格轉化得到的有限元網(wǎng)格為六節(jié)點五面體的退化六面體網(wǎng)格(見圖6)。

圖6 三角形內(nèi)表面網(wǎng)格拉伸圖

以圖6為例，推導同一節(jié)點連接的所有單元的合矢量的生成公式及拉伸距離：

(4)

(5)

(6)

(7)

(8)

(9)

根據(jù)矢量歸一化特性有：

(10)

(11)

式(11)中：D表示物體厚度；DL表示合矢量負方向的截距。

(12)

由式(10)—式(12)得：

(13)

由式(10)—式(13)可求得內(nèi)表面網(wǎng)格點的坐標：

(14)

式(14)中：x、y、z分別表示外表面三角網(wǎng)格點坐標；x′、y′、z′分別表示內(nèi)表面三角網(wǎng)格點的坐標。

2 變尺寸等厚度STL表面模型的3D打印數(shù)據(jù)處理方法

通過理論研究得到單一厚度STL表面模型生成有限元網(wǎng)格的方法，在此基礎上提出變尺寸等厚度STL表面模型的3D打印數(shù)據(jù)處理方法。

變尺寸等厚度結構三維模型如圖7所示。在處理STL文件得到滿足轉化條件的外表面網(wǎng)格后，沿著合矢量負方向進行拉伸，不同于單一厚度結構外表面拉伸，其外表面要先沿著合矢量負方向拉伸一個厚度hL1，沿合矢量負方向得到新節(jié)點p1，再沿著合矢量負方向拉伸一個距離hL2-hL1，沿合矢量負方向得到新節(jié)點p2，過新節(jié)點p2向邊L作垂線生成垂足點p3，點p3取代點p2與其他節(jié)點相連接，生成變尺寸等厚度網(wǎng)格圖如圖8所示。

圖8 變尺寸等厚度網(wǎng)格圖

求變尺寸等厚度結構同一節(jié)點連接的所有單元的合矢量及拉伸距離公式：

(15)

(16)

根據(jù)適量歸一化特性有：

(17)

(18)

需拉伸的距離：

(19)

式(19)中：h1表示為圖8中左側較薄部分的物體厚度；hL1表示為合矢量負方向的截距，即生成退化六面體網(wǎng)格需沿合矢量負方向拉伸的距離。

(20)

式(20)中：h2表示為圖8中右側較厚部分的物體厚度；hL2表示為合矢量負方向的截距；hL2-hL1為表面需要拉伸的距離。

由式(20)可求得變尺寸等厚度結構內(nèi)表面網(wǎng)格點p1、p2、p3的坐標。p1點坐標：

(21)

式(21)中：x、y、z分別表示變尺寸等厚度結構外表面三角網(wǎng)格點坐標；x′1、y′1、z′1分別表示變尺寸等厚度結構內(nèi)表面p1點的坐標。

p2點坐標：

(22)

式(22)中：x、y、z分別表示變尺寸等厚度結構外表面三角網(wǎng)格點坐標；x′2、y′2、z′2分別表示變尺寸等厚度結構內(nèi)表面p2點的坐標。

p3點坐標：

(23)

式(23)中：x、y、z分別表示變尺寸等厚度結構外表面三角網(wǎng)格點坐標；x′3、y′3、z′3分別表示變尺寸等厚度結構內(nèi)表面p3點的坐標。

3 實例計算

為驗證表面類STL網(wǎng)格轉化為有限元計算網(wǎng)格的適用性，以某水電站彎肘式尾水管為例，進行了如下計算嘗試。實例計算選用的水電站彎肘式尾水管，其直錐段上部半徑1.65 m,下部半徑2.33 m，高度3.66 m；彎肘段上部為圓結束處為方，上部中心到彎管結束水平距離6 m，高4.46 m；擴散段為入口到出口水平距離為8.85 m。入口出口均為矩形，入口長8.98 m，寬2.89 m；出口長8.98 m，寬4.46 m。

3.1 尾水管模型建立及有限元網(wǎng)格生成

得到3D打印格式的STL文件，通過程序轉化為有限元網(wǎng)格，最后導入到有限元軟件中計算。通過程序生成得到的有限元法表示的部分點坐標和單元信息的文件如圖9所示，基于該文件可生成不同格式的有限元文件，不同格式的有限元文件可導入不同的有限元軟件中計算，本文將有限元文件導入ANSYS中，生成的有限元網(wǎng)格如圖10所示。

圖9 文件中部分節(jié)點

圖10 尾水管有限元網(wǎng)格細部圖

3.2 水電站尾水管有限元計算

在ANSYS Workbench中進行有限元計算，力學參數(shù)及工況設定如下：

(1)設定模型邊界條件：尾水管底部固定。

(2)材料力學參數(shù)：尾水管材料為鋼；密度為7 890 kg/m3；彈性模量為2.06×1011Pa；泊松比為0.27。

(3)在工況1情況下，對尾水管施加一個自重荷載；工況2：對尾水管同時施加一個自重荷載和靜水壓力。工況1下計算得到尾水管的總變形圖、正應力圖、等效應力圖見圖11-圖13；在工況2情況下，計算得到總變形圖見圖14-圖16。由圖11-圖16可知，尾水管的應力及變形情況均在正常情況范圍內(nèi)。

圖11 尾水管總變形圖

圖12 尾水管Z向應力圖

圖13 尾水管等效應力圖

圖14 尾水管總變形圖

圖15 尾水管Z向應力圖

圖16 尾水管等效應力圖

4 結 語

本文基于激光掃描或攝影成像獲取的表面模型進行3D打印數(shù)據(jù)格式處理及有限元網(wǎng)格的轉化，在單層表面STL網(wǎng)格的基礎上得到單一厚度的3D打印數(shù)據(jù)格式，并提出了可變厚度的STL表面模型3D打印數(shù)據(jù)及有限元網(wǎng)格的新方法。結合某水電站尾水管建立單層表面模型進行有限元網(wǎng)格轉化，設定尾水管受自重和靜水壓力兩個工況下進行有限元計算，得到尾水管變形、應力狀態(tài)。通過這一尾水管實例，對于復雜水工結構利用表面類STL網(wǎng)格轉化為有限元網(wǎng)格，探索出一種新的途徑。