基于自適應三角化的船舶曲面分段網格細分

李純金，楊秋林

(江蘇科技大學 機械工程學院，江蘇 鎮江 212003)

0 引 言

船舶外板曲面展開一般通過三角形網格或者四邊形網格進行展開。曲面網格劃分成為外板曲面展開的重要部分，對船舶外板曲面展開的精度有著重要的影響。網格過于稀疏會使展開誤差加大，過密則會降低展開效率。對于給定的船舶型值表數據，重構外板曲面后，曲面網格的劃分引起了廣泛的關注。在世博軸“陽光谷”工程設計中，李承銘[1]先對參數域曲面進行人工細分后得到平面網格，之后將二維平面網格映射回三維空間曲面，實現了對該參數曲面的網格劃分，網格劃分效果較好，但其步驟繁瑣，通用性不高。丁慧[2]基于傳統的映射方法，基于等參數線分割法的劃分方法生成自由曲面網格，通過等弦長分割方式對曲面參數化不均勻的問題提出初步解決方案，但其生成的網格曲面單元仍存在劃分混亂情況。江存[3]提出自定義單元法的自由曲面網格劃分，建立網格拓撲的參數化表達，改善了網格映射時的網格失真，但其方法試用范圍單一，對非NURBS曲面不適用。王結臣[4]通過柵格對空間散亂點集進行劃分，以柵格為整體進行區域網格劃分。陳永就[5]與衛洪春[6]分別通過直接剖分及遞歸程序分別對四邊形及三角形網格進行了處理，網格擴展過程中變形量較小。石磊[7]對現有的細分算法做出了分析，并同時指出了現有細分算法的劣勢。針對船舶外板曲面的網格細分算法研究不多，由于船舶外板曲面曲率比較復雜，需考慮網格劃分疏密對于以后外板曲面展開精度的影響。

1 Delaunay 三角網格化

1.1 Delaunay定義

Delaunay三角網格劃分為數值分析及圖形學方面重要技術手段，如圖1所示，Delaunay三角網格是將空間測量數據點投影到平面來實現的二維劃分法。算法基本思想為每一個Delaunay三角形的外接圓內部不包含其他結點，一旦檢測到三角形外接圓含有其他結點后，就必須修改局部的網格劃分，直至滿足起始條件為止。

1.2 Delaunay三角網格劃分算法

Delaunay三角網格劃分步驟為：

1）構造大外接圓，使其包含所有結點，進行步驟2；

2）每次導入一個結點，重復進行步驟3～步驟5；

3）找出已有三角形中哪些外接圓包含新引入的結點，轉步驟4；

4）刪除這些三角形中離新結點最近的一條邊，轉步驟5；

5）將新結點與老結點連接，產生新的三角形，算法示意圖如圖2所示。

多邊形域由Delaunay三角網格劃分后生成的網格還不能滿足幾何造型或者有限元分析的需求，為此需對初步劃分后的網格作進一步的細化，在邊界或者圖形區域內部插入適當的節點，提高網格的質量。

2 細分曲面算法

對于給定的點云數據或者三角網格，構造細分曲面使得曲面插值于點云或者三角網格頂點，就是細分曲面插值問題。

2.1 蝶形算法

蝶形算法為一種插值細分算法，如圖3所示在每一條邊上引入一個新頂點，頂點的值由模型中8個頂點的值加權平均得到。權重系數設置為：

其中為張力系數 ，與表面的光滑程度有關。蝶形算法產生C1連續的細分網格，但是對奇異點的處理不夠。

2.2 多面體算法

多面體算法在三角形網格細分中屬于最簡便的算法，算法基本思想是在三角形網格每2個頂點連線中間引入新的頂點，將原有單個三角形網格細分為4個三角網格。新引入頂點的位置通過取其所在連線兩端2個頂點加權平均值確定，該算法不受奇異點影響，常用于細分算法的原理介紹。圖4為多面體細分示意圖。

2.3 Loop算法

Loop算法是一種逼近細分算法，其算法步驟主要分為3步，首先生成邊界點，之后移動原有頂點，最后生成光滑表面（切平面連續）。

若以為曲面中心頂點，為相鄰頂點構成傘狀網格，再根據Loop算法就可得到下式：

綜合上述兩式得到如下局部細分矩陣形式：

中心頂點處的局部細分矩陣即為式（4）中的矩陣，采用靜態的網格細分模式，與網格頂點在網格中所處位置無關，只與階數有關。

3 自適應三角網格劃分算法（ATM）

3.1 ATM算法思想

自適應網格劃分（Adaptive Triangular Meshing）算法通過自動估算曲面曲率的大小，結合Delaunay算法，自動選擇合適的劃分密度，避免了人工劃分的繁瑣，提高了劃分效率。針對船舶外板曲面自適應網格劃分算法基本步驟如下：

1）提取船舶型值表數據，通過B樣條工具進行參數化曲面重構，得到參數化曲面；

2）選定曲面適當區域，在參數化曲面方向上選取相向的2條曲線，從一端的頂點開始搜索，計算步長兩端點弦長或者弧長（估算值）差的絕對值，將其值稱為弦弧差；

3）設定弦弧差許用范圍，將步驟2得到的值與弦弧差范圍比較，如果在弦弧差范圍之內，記錄步長，并生成劃分用網格線，否則步驟4；

4）調整步長重新判斷，直至滿足弦弧差范圍；

5）重復步驟2～步驟4，直到搜索至邊界的終點；

6）所有搜索完畢后得到的是四邊形網格，此時進行Delaunay三角化生成細分網格。

3.2 步長劃分準則

3.1節所述ATM算法步驟適用于矩形參數域。對于非矩形參數域而言，需要先對曲面進行擴展補充，通過方 向的最大與最小參數值確定一個矩形域，此矩形域為原先非矩形參數域的最小外接矩形域。

以向網格劃分為例，如圖5所示，原先參數域為非矩形參數域，建立如圖的坐標系，通過最大參數值及最小參數值構建最小外接矩形參數域，包含原有的多邊形域中所有點。B與b之間距離的絕對值為劃分的步長值，將其與設定的弦弧差范圍比較。

選中由最值形成的參數化區域，由a、b兩點開始搜索，給定初始步長l，求出a與A，b和B之間的弦弧差的值，其中較大的弦弧差的值記作，較小的記作。 給定的許用弦弧差最小值為，許用弦弧差最大值為。判斷當則滿足條件，步長合適并生成劃分網格線，否則需要重新修改步長再判定。

3.3 步長關系調整

步長關系調整基于精度為主，效率為輔的原則。實際與許用弦弧差關系總共分為7種，如圖6所示，上面部分圖形為實際弦弧差絕對值，左端為最小值，右端為最大值，下面部分圖形為許用弦長差絕對值，左端為許用弦弧差范圍最小值，右端為許用弦弧差最大值。

針對7種步長關系圖進行如下調整：

1）在情況1中，實際的弦弧差范圍均處于許用弦弧差范圍內，，，滿足判斷條件，可以記錄步長，生成網格劃分線，并進行下一步搜索；

2）在情況2中，實際弦弧差與許用弦弧差關系為，，此時需要加大步長，在其滿足判斷條件后記錄步長，生成網格劃分線；

3）在情況3中，實際弦弧差最大值小于許用弦弧差最大值，關系為，，跟情況2的不同之處在于此時實際弦弧差最大值接近于許用弦弧差，不能加大步長，記錄步長，生成網格劃分線；

4）情況4中實際弦弧差與許用弦弧差關系為，，此時需要縮小步長，使得，即滿足情況1或情況3的狀態，當滿足這個條件時記錄步長，并生成網格劃分線；

5）情況5中實際弦弧差與許用弦弧差關系為，，此時需要增大步長值，在其滿足判斷條件后記錄步長，生成網格劃分線；

6）情況6中實際弦弧差與許用弦弧差關系為，此時需要縮小步長，記錄步長，在其滿足判斷條件后記錄步長，生成網格劃分線；

7）情況7中實際弦弧差與許用弦弧差關系為，，此時減少步長值，使其達到情況3的狀態。

4 實例應用

算法實現電腦配置為PC i3CPU，M350 @2.27 GHz RATM4G。根據提供的船舶型值表數據首先通過B樣條進行曲面重構得到其參數化曲面，之后通過自適應算法得出船舶曲面外板實現網格劃分效果。具體實現流程如圖7所示。

船舶外板曲面實例效果如圖8所示。經過自動調整劃分的網格對于船舶外板的復雜曲面展開的精度更有益處。

下式為計算曲面展開面積誤差的表達式：

其中和為展開前后的三角網格的對應三角形的面積值；為空間曲面第個頂點相連的三角形的個數。

如表1所示，在對不同劃分網格數的曲面進行展開后，得到的曲面面積變化誤差也不同。對初始曲面分別進行劃分處理，得到網格數為460，900，1 840的三角形網格曲面以及自適應細分后的1 386個三角片的曲面，借助現有的三角形展開方法進行展開后，得到相應的面積誤差。可見，對于精細劃分后的曲面展開來講，其展開的面積誤差較小，對于船板制造加工具有指導意義。

表1 不同網格數曲面展開誤差對比Tab.1 The number of different grid surface expansion error comparison

5 結 語

本文取船舶曲面分段進行曲面網格的劃分研究，針對船舶曲面外板曲面的復雜特征而言，選擇合適的網格劃分方式對外板展開尤為重要。通過B樣條工具進行曲面重構，之后對重構出的四邊形網格作進一步的細分網格處理，并通過曲面展開驗證細分網格的效果?，F有的網格細分方法對劃分網格的準確性存在一定缺陷，將Delaunay三角網格劃分與自適應劃分算法結合起來，利用自適應劃分網格的準確性，實現船舶外板曲面的三角網格化處理，為曲面展開做好鋪墊。

[1]李承銘, 盧旦.自由曲面單層網格的智能布局設計研究[J].土木工程學報, 2011, 3: 1–7.

[2]丁慧.自由形態空間網格結構的網格設計方法研究與實現[D].杭州: 浙江大學, 2014.

[3]江存.自由曲面空間網格結構網格劃分、優化及力性能研究[D].杭州: 浙江大學, 2015.

[4]陳永就, 林川.基于直接剖分的球面四邊形離散格網生成方法[J].地理空間信息, 2015, 3: 130–132+12.

[5]衛洪春.多邊形三角剖分與三角細分的研究與實現[J].計算機與現代化, 2015, 7: 65–68, 76.

[6]石磊, 薛珊.基于三角形網格的幾種典型細分曲面方法概述[J].赤峰學院學報(自然科學版), 2013, 14: 13–14.

[7]LI Qi-rui.Surfaces expanding by the power of the Gausscurvature flow[J].Proceedings of the American Mathematical Society, 2010.