陽(yáng)湘安，阮 鋒

（1.廣東技術(shù)師范學(xué)院 機(jī)電學(xué)院，廣東 廣州510635；2.華南理工大學(xué) 機(jī)械與汽車工程學(xué)院，廣東 廣州510640）

0 引 言

在薄板的有限元成形過(guò)程仿真中，經(jīng)常采用四邊形網(wǎng)格或三角形網(wǎng)格或兩者的混合網(wǎng)格對(duì)模型進(jìn)行曲面的離散化。其中，離散網(wǎng)格頂點(diǎn)法矢對(duì)于計(jì)算頂點(diǎn)曲率、網(wǎng)格光順、網(wǎng)格優(yōu)化和網(wǎng)格變形以及離散網(wǎng)格的曲面擬合等都具有重要的作用［1－3］，對(duì)于有限元的計(jì)算結(jié)果也有一定的影響。

計(jì)算離散網(wǎng)格頂點(diǎn)法矢的主要方法目前可以區(qū)分為2類，即基于二次誤差矩陣法向表決的法矢計(jì)算方法和權(quán)重平均的法矢計(jì)算方法。第1類計(jì)算方法精度較高，但是計(jì)算的法矢仍然依賴于某權(quán)重，同時(shí)，對(duì)于兩面偏差角大于或等于90°時(shí)，需要調(diào)整算法，甚至得不到結(jié)果，而且通過(guò)測(cè)地線來(lái)計(jì)算頂點(diǎn)的鄰域也很耗時(shí)［4，5］。對(duì)于曲面的離散網(wǎng)格頂點(diǎn)的法矢計(jì)算，更多的是采用第2類計(jì)算方法［6－8］，考慮到網(wǎng)格的不規(guī)則性，也有提出計(jì)及質(zhì)心位置坐標(biāo)的權(quán)重方法或者單獨(dú)引入網(wǎng)格的形狀因子來(lái)計(jì)算權(quán)重［9，10］，但是，對(duì)于曲面的有限元網(wǎng)格劃分，由于有限元計(jì)算的要求，曲面劃分后的離散網(wǎng)格的形狀通常是較為規(guī)整的，因此，需要正確評(píng)估這種情況下各種權(quán)重平均的法矢計(jì)算方法的計(jì)算精度和效率，同時(shí)，目前針對(duì)三角形網(wǎng)格的頂點(diǎn)法矢算法比較成熟，但關(guān)于四邊形網(wǎng)格、四邊形網(wǎng)格與三角形網(wǎng)格并存的混合網(wǎng)格的頂點(diǎn)法矢計(jì)算的研究則較少，基于這種情況本文提出了解決這一問(wèn)題的算法。

1 三角形離散網(wǎng)格曲面的法矢計(jì)算

1.1 基于二次誤差矩陣法向表決的法矢計(jì)算方法

這種方法的思想首先是確定頂點(diǎn)的鄰域，即一系列的三角形網(wǎng)格平面fi，則點(diǎn)x 到三角形網(wǎng)格平面的距離為

式中：p——三角形網(wǎng)格平面上的某個(gè)點(diǎn)，δi＝－為三角形網(wǎng)格平面到原點(diǎn)之間的距離，——三角形網(wǎng)格平面的法矢。

而帶權(quán)重wi的n 個(gè)通過(guò)同一個(gè)點(diǎn)p 的平面與x 的距離的平方為

這里采用權(quán)重wi是考慮各三角形網(wǎng)格平面的貢獻(xiàn)差異，而采用距離的平方則考慮到對(duì)方程求最小值后得到二次誤差矩陣方程，對(duì)于所求網(wǎng)格頂點(diǎn)，其與1－ring三角形網(wǎng)格平面的距離為0。式 （2）最小化后得到如下方程

一般來(lái)說(shuō)，A 是對(duì)稱、半正定的，假設(shè)求得其特征值為λi（λ1≥λ2≥λ3），對(duì)應(yīng)的特征向量為，則其中的特征向量即為要求取的頂點(diǎn)法矢。

1.2 權(quán)重平均的法矢計(jì)算方法

基于二次誤差矩陣法向表決的法矢計(jì)算方法需要較大的工程計(jì)算量，相比來(lái)說(shuō)，權(quán)重平均的法矢計(jì)算方法則簡(jiǎn)便得多，然而，很多確定權(quán)重的方法都是基于經(jīng)驗(yàn)來(lái)進(jìn)行的，只對(duì)某類曲面的網(wǎng)格有較高精度。典型的權(quán)重確定方法有：

（1）等權(quán)重方法，即各三角形網(wǎng)格的權(quán)重相等且一般都取為1，因而，網(wǎng)格頂點(diǎn)法矢的計(jì)算為

（2）面積權(quán)重法。這種方法認(rèn)為三角形網(wǎng)格的面積對(duì)頂點(diǎn)法矢的影響是主要的，因而以各三角形網(wǎng)格的面積Ai作為各自法矢的權(quán)重，網(wǎng)格頂點(diǎn)法矢的計(jì)算則為

（3）角度權(quán)重法。這種方法則認(rèn)為頂點(diǎn)相鄰兩邊的夾角對(duì)頂點(diǎn)法矢的影響具有最大影響，因而根據(jù)頂點(diǎn)相鄰兩邊的夾角αi大小來(lái)確定權(quán)重，網(wǎng)格頂點(diǎn)法矢的計(jì)算為

（4）角度和面積乘積的權(quán)重法。在這一方法中，權(quán)重的確定綜合考慮了角度和面積的影響，其表達(dá)式為

（5）質(zhì)心位置的權(quán)重確定方法。這種方法認(rèn)為質(zhì)心的位置能夠比較全面地反映三角形網(wǎng)格的形狀。為了考慮三角形網(wǎng)格形狀對(duì)權(quán)重確定的影響，權(quán)重的確定采用如下表達(dá)式

式中：ωi——質(zhì)心權(quán)重，gi——三角形網(wǎng)格的質(zhì)心位置。

1.3 Max方法的原理

1.2節(jié)中各種權(quán)重平均的法矢計(jì)算方法都是根據(jù)經(jīng)驗(yàn)來(lái)確定的，缺乏相應(yīng)的理論依據(jù)，而Max算法雖然形式上體現(xiàn)為一種通過(guò)權(quán)重平均來(lái)確定網(wǎng)格的頂點(diǎn)法矢，但本質(zhì)上來(lái)說(shuō)是一種通過(guò)外接球面擬合四面體來(lái)確定頂點(diǎn)法矢的方法，因而具有相應(yīng)的理論依據(jù)。

圖1所示，對(duì)于以Q 為頂點(diǎn)、Vi為節(jié)點(diǎn)的三角形網(wǎng)格系列，選取其中的2 個(gè)三角形網(wǎng)格QV1V2和QV2V3，則點(diǎn)Q、V1、V2和V3構(gòu)成一個(gè)四面體，而通過(guò)四面體的4個(gè)頂點(diǎn)可以做出唯一的一個(gè)外接球面，于是，通過(guò)三角形網(wǎng)格邊線QV1、QV2、QV3的中點(diǎn)A、B、C可以構(gòu)建其相應(yīng)的3個(gè)法平面方程，對(duì)3個(gè)方程的聯(lián)立求解則可以得到球面圓心O 的坐標(biāo)，而矢量O→Q 則為由這2個(gè)相鄰三角形網(wǎng)格確定的法矢。

圖1 Max算法原理

對(duì)網(wǎng)格頂點(diǎn)及其鄰接n個(gè)三角形網(wǎng)格循環(huán)采用上述方法并將結(jié)果歸一化，則得到最終的網(wǎng)格頂點(diǎn)法矢

2 有限元網(wǎng)格頂點(diǎn)的法矢算法

2.1 典型曲面樣本的有限元網(wǎng)格模型

為了比較各種權(quán)重平均的法矢計(jì)算方法的結(jié)果精度，選取以下4種曲面為對(duì)象，通過(guò)有限元網(wǎng)格劃分后得到四邊形網(wǎng)格模型或者同時(shí)帶有四邊形網(wǎng)格和三角形網(wǎng)格的混合網(wǎng)格模型。劃分后的網(wǎng)格模型如圖2所示。

圖2 典型曲面樣本有限元網(wǎng)格模型

在上述曲面的有限元網(wǎng)格模型中，半球面和圓錐面大部分為四邊形網(wǎng)格，但在邊界處出現(xiàn)三角形網(wǎng)格，這種情況在實(shí)際中最為常見(jiàn)，而環(huán)面和圓柱面則由于比較規(guī)整，因而劃分出的網(wǎng)格只有四邊形網(wǎng)格。由于實(shí)際中三角形網(wǎng)格的頂點(diǎn)法矢計(jì)算較為成熟，因而，可以采取將四邊形網(wǎng)格轉(zhuǎn)化為三角形網(wǎng)格的頂點(diǎn)法矢計(jì)算模式。

2.2 有限元網(wǎng)格頂點(diǎn)法矢的算法流程

將四邊形網(wǎng)格轉(zhuǎn)化為三角形網(wǎng)格的方法有很多種，典型的有3種，如圖3所示。方案 （a）中，對(duì)于四邊形網(wǎng)格中的某節(jié)點(diǎn)Vij，將Vij與Vi－1，j－1和Vi＋1，j＋1連接 起來(lái) 即 產(chǎn)生新的邊，這樣四邊形網(wǎng)格就分解成了三角形網(wǎng)格；方案（b）中，對(duì)四邊形網(wǎng)格中的滿足i＋j為偶數(shù)的點(diǎn)Vij，將點(diǎn)Vij與其1－鄰域中的對(duì)角點(diǎn)連接，產(chǎn)生的線段作為加入的新邊，這樣也形成了三角形網(wǎng)格；方案 （c）中，則直接對(duì)矩形域上的所有頂點(diǎn)作Delaunay三角剖分［11］。

圖3 四邊形網(wǎng)格到三角形網(wǎng)格轉(zhuǎn)換的典型方法

對(duì)于網(wǎng)格頂點(diǎn)法矢的計(jì)算來(lái)說(shuō)，應(yīng)該盡量利用頂點(diǎn)周圍鄰近節(jié)點(diǎn)的信息，而上述方案都沒(méi)有滿足這一要求，并且，由于物理上并不需要將四邊形網(wǎng)格劃分為三角形網(wǎng)格，因而，對(duì)于所需求取法矢的網(wǎng)格頂點(diǎn)，可以首先搜尋到與其直接相鄰的其他網(wǎng)格頂點(diǎn)，將其中屬于同一個(gè)四邊形網(wǎng)格單元的網(wǎng)格頂點(diǎn)連接起來(lái)，從而與所需求取法矢的網(wǎng)格頂點(diǎn)一起構(gòu)成計(jì)算所需的三角形網(wǎng)格。如對(duì)于頂點(diǎn)Vij，將與其直接相鄰的節(jié)點(diǎn)Vi，j＋1、Vi－1，j、Vi，j－1和Vi＋1，j依次連接起來(lái)構(gòu)成四條新邊，進(jìn)而與頂點(diǎn)Vij共同構(gòu)成計(jì)算其頂點(diǎn)法矢的三角形網(wǎng)格，如圖4所示。

圖4 四邊形網(wǎng)格到三角形網(wǎng)格的處理

編程計(jì)算流程為首先將四邊形網(wǎng)格轉(zhuǎn)換為三角形網(wǎng)格模式，再針對(duì)三角形網(wǎng)格采用不同的權(quán)重平均算法來(lái)計(jì)算，具體流程如圖5所示。

圖5 有限元網(wǎng)格頂點(diǎn)法矢的算法流程

采用各種權(quán)重平均的法矢計(jì)算方法對(duì)上述曲面樣本網(wǎng)格模型分別進(jìn)行編程，計(jì)算各頂點(diǎn)的法矢，得到的結(jié)果見(jiàn)表1。

各類型曲面網(wǎng)格頂點(diǎn)法矢計(jì)算結(jié)果的最小均方根偏差和最小絕對(duì)值平均偏差對(duì)應(yīng)的計(jì)算方法見(jiàn)表2。

3 結(jié)果分析

從表1和表2可以看出，對(duì)于半球面、環(huán)面和圓柱面，Max方法相對(duì)于其他權(quán)重平均方法來(lái)說(shuō)，其計(jì)算結(jié)果的均方根偏差和絕對(duì)值平均偏差都最小，對(duì)于圓錐面來(lái)說(shuō)，其計(jì)算結(jié)果的絕對(duì)值平均偏差也最小，而均方根偏差為0.0046 mm，相對(duì)于其他權(quán)重平均方法的最小值0.0042 mm 來(lái)說(shuō)，差距也并不是很大。

這說(shuō)明，Max 方法不但計(jì)算精度很高，而且具有很強(qiáng)的適應(yīng)性。這主要是因?yàn)镸ax方法在計(jì)算頂點(diǎn)法矢時(shí)對(duì)于兩相鄰的三角形網(wǎng)格構(gòu)成的四面體采用外接球面進(jìn)行擬合，精確計(jì)算出2個(gè)三角形網(wǎng)格時(shí)頂點(diǎn)的法矢，并循環(huán)采用這一方法計(jì)算出其它兩相鄰三角形網(wǎng)格時(shí)的頂點(diǎn)法矢，由于這種方法在計(jì)算時(shí)同時(shí)考慮了兩個(gè)相鄰三角形網(wǎng)格的法矢，因此，當(dāng)將所有結(jié)果作矢量加得到最終頂點(diǎn)的法矢后，其結(jié)果精度必然就會(huì)很高；而且，由于對(duì)兩相鄰的三角形網(wǎng)格構(gòu)成的四面體采用外接球面進(jìn)行擬合，因而也必然計(jì)及了三角形網(wǎng)格的形狀規(guī)整性，這一結(jié)論從文獻(xiàn)［10］中也可以看出。從計(jì)算效率來(lái)說(shuō)，由于Max方法最終可以推導(dǎo)為如式 （9）所示的簡(jiǎn)單形式，因而其算法的效率也不會(huì)很低。相對(duì)而言，其它的權(quán)重方法只是考慮如何對(duì)一個(gè)單獨(dú)的三角形網(wǎng)格確定權(quán)重，而且其權(quán)重的確定也主要是從經(jīng)驗(yàn)的角度來(lái)考慮，沒(méi)有相應(yīng)的理論依據(jù)，因而其結(jié)果精度不但較低，而且對(duì)于不同類型的曲面其結(jié)果精度不穩(wěn)定。

表1 各方法的偏差比較

表2 最小偏差對(duì)應(yīng)的計(jì)算方法

從表1和表2還可以看到，采用等權(quán)重法計(jì)算后結(jié)果的均方根偏差和絕對(duì)值平均偏差也都比較小，其中對(duì)于圓錐面來(lái)說(shuō)，其均方根偏差和絕對(duì)值平均偏差都最小，對(duì)于半球面，其均方根偏差也最小。究其原因，主要是由于有限元計(jì)算的要求，通常劃分后的網(wǎng)格都比較規(guī)整，而其它的經(jīng)驗(yàn)性權(quán)重平均方法中引入的權(quán)重因子由于并不能全面反映三角形網(wǎng)格的形狀差異，因而這些方法并不能有效提高計(jì)算結(jié)果的精度。由于等權(quán)重的計(jì)算簡(jiǎn)單，又具備一定的計(jì)算結(jié)果精度，因而，對(duì)于有限元網(wǎng)格頂點(diǎn)的法矢計(jì)算要求不是非常高的情況下，從計(jì)算效率的角度考慮是可以采用的。

4 結(jié)束語(yǔ)

文中首先闡述了計(jì)算三角形網(wǎng)格頂點(diǎn)法矢的各種方法的原理和特點(diǎn)，針對(duì)典型曲面的有限元網(wǎng)格的頂點(diǎn)法矢計(jì)算，分析了基于經(jīng)驗(yàn)的權(quán)重平均法與Max方法的計(jì)算結(jié)果精度。結(jié)果表明，對(duì)于計(jì)算結(jié)果精度要求較高且計(jì)算效率要求也較高的情況，可以采用等權(quán)重法來(lái)確定曲面有限元網(wǎng)格的頂點(diǎn)法矢，而對(duì)于計(jì)算結(jié)果精度要求很高的情況，則應(yīng)該采用基于外接球面擬合四面體的Max 方法達(dá)到這一目標(biāo)。

文中針對(duì)四邊形網(wǎng)格以及四邊形和三角形網(wǎng)格并存的情況，提出了將四邊形網(wǎng)格處理為三角形網(wǎng)格并有效計(jì)算出網(wǎng)格頂點(diǎn)法矢的方法，這一方法盡量考慮到了利用網(wǎng)格頂點(diǎn)相鄰節(jié)點(diǎn)的的信息；設(shè)計(jì)的算法不但可以單獨(dú)處理四邊形網(wǎng)格或三角形網(wǎng)格，而且能夠處理曲面有限元網(wǎng)格劃分后四邊形網(wǎng)格和三角形網(wǎng)格混合并存的情況。

［1］HE Qiang，ZHANG Shusheng，BAI Xiaoliang.Mesh smoothing algorithm based on local surface approximation ［J］.Journal of Harbin Institute of Technology，2011，43 （5）：89－93（in Chinese）.［賀強(qiáng)，張樹(shù)生，白曉亮.基于局部曲面逼近的網(wǎng)格光順?biāo)惴?［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2011，43 （5）：89－93.］

［2］CHEN Zhiyang，HAN Chunlei.Simplification and reconstruction of mesh model in CAE ［J］.Computer Systems ＆Applications，2011，20 （5）：188－191 （in Chinese）. ［陳志楊，韓春雷.面向CAE的網(wǎng)格模型特征簡(jiǎn)化與重建 ［J］.計(jì)算機(jī)系統(tǒng)應(yīng)用，2011，20 （5）：188－191.］

［3］ZHANG Aiwu.Surface reconstruction with normal control［D］.Jinan：Shandong University，2009 （in Chinese）.［張愛(ài)武.法矢控制的網(wǎng)格曲面重建若干問(wèn)題研究 ［D］.濟(jì)南：山東大學(xué)，2009.］

［4］Di Angelo L，Di Stefano P.Experimental comparison of methods for differential geometric properties evaluation in triangular meshes［J］.Computer－Aided Design and Applications，2011，8 （2）：193－210.

［5］Hyun Soo Kim，Han Kyun Choi，Kwan H Lee.Feature detection of triangular meshes based on tensor voting theory ［J］.Computer－Aided Design，2009，41 （1）：47－58.

［6］YANG Nan，XIAO Jiangchao，WANG Minghai.Estimation of normal and curvature based on the triangle mesh model［J］.Modern Manufacturing Engineering，2010，38 （3）：104－107 （in Chinese）.［楊楠，校江超，王明海.基于三角網(wǎng)格模型的法矢及曲率估算［J］.現(xiàn)代制造工程，2010，38 （3）：104－107.］

［7］Kim Hyoung－Seok，Kim Ho－Sook.New computation of normal vector and curvature ［J］.Wseas Transactions on Computers，2009，8 （10）：1661－1670.

［8］ZHANG Yajuan，LI Ze’an，CHENG Chen.Modified smoot－h(huán)ing algorithm of triangular mesh model［J］.Computer Engineering，2012，38 （9）：220－222 （in Chinese）. ［章雅娟，李澤安，程晨.一種改進(jìn)的三角網(wǎng)格模型光順?biāo)惴?［J］.計(jì)算機(jī)工程，2012，38 （9）：220－222.］

［9］Chen Yueping，Gao Jian，Wen Hao，et al.Estimation normal vector of triangular mesh vertex by angle and centroid weights and its application ［J］.TELKOMNIKA，2013，11 （4）：1841－1848.

［10］PENG Yuhui，GAO Chenghui.An improved algorithm for vertex normal computation of triangular meshes based on shape correction ［J］.Journal of Image and Graphics，2010，15（1）：142－148 （in Chinese）.［彭育輝，高誠(chéng)輝.基于形狀修正的三角網(wǎng)格模型頂點(diǎn)法矢估算方法 ［J］.中國(guó)圖象圖形學(xué)報(bào)，2010，15 （1）：142－148.］

［11］WANG Cheng’en，CUI Dongliang，YAN Zhiyuan，et al.Finite element triangle mesh generation in planar area ［J］.Computer Integrated Manufacturing Systems，2011，17 （2）：256－260 （in Chinese）.［王成恩，崔東亮，嚴(yán)智遠(yuǎn)，等.平面區(qū)域有限元三角形網(wǎng)格劃分研究及實(shí)現(xiàn) ［J］.計(jì)算機(jī)集成制造系統(tǒng)，2011，17 （2）：256－260.］