郭新強， 張見明， 覃先云， 穆東輝， 宋 喬

（1. 湖南大學(xué)汽車車身先進(jìn)設(shè)計制造國家重點實驗室，湖南 長沙 410082；2. 北京第二機床廠有限公司，北京 100072）

曲面網(wǎng)格生成的方法通常分為直接法和映射法[1]。直接法是指曲面網(wǎng)格的劃分直接在曲面的物理空間進(jìn)行；映射法是首先利用平面域網(wǎng)格生成方法在曲面的二維參數(shù)空間進(jìn)行網(wǎng)格剖分，然后將剖分結(jié)果映射到三維物理空間形成曲面網(wǎng)格。本文利用反映曲面扭曲程度的黎曼度量[2-4]，將用于生成平面四邊形網(wǎng)格技術(shù)的鋪磚法[5]擴展運用于任意復(fù)雜曲面的網(wǎng)格生成。鋪磚法首先由Blacker T D和 Stephenson M B提出的，該法具有較高的邊界靈敏度，可以較好的擬合網(wǎng)格區(qū)域復(fù)雜的邊界形狀，能生成質(zhì)量較高的網(wǎng)格。在原有的鋪磚法中，首先沿著網(wǎng)格邊界生成一排單元，然后再考慮相交情況，使得相交處理復(fù)雜。在本文實現(xiàn)該方法過程中，每生成一個單元時就預(yù)先判斷相交情況，根據(jù)判斷對相交單元做簡單的處理，省去了原方法中對多次相交優(yōu)先處理環(huán)節(jié)，從而使得整個的相交處理變得簡單。

本文利用改進(jìn)的鋪磚法生成的四邊形單元是用于一種新的結(jié)構(gòu)分析算法——邊界面法[6]前處理的二維參數(shù)曲面網(wǎng)格。邊界面法由張見明等人最近提出，該法直接利用CAD實體模型的邊界參數(shù)曲面信息在二維參數(shù)空間中劃分網(wǎng)格并在曲面的二維參數(shù)空間內(nèi)進(jìn)行邊界積分和變量插值。在數(shù)值積分過程中，積分點的幾何數(shù)據(jù)直接由參數(shù)曲面的參數(shù)變量計算得到[7]，而不是基于單元由多項式插值近似的，避免了幾何誤差。該方法中所使用的四邊形單元定義在曲面的二維參數(shù)空間，而不是像有限元法中的單元定義在三維物理空間。有限元法的幾何模型是通過單元插值近似的，當(dāng)網(wǎng)格較稀疏時具有明顯的幾何誤差。邊界面法是利用二維參數(shù)空間的值進(jìn)行數(shù)值計算，現(xiàn)有的有限元前處理商業(yè)軟件能輸出模型的三維物理空間坐標(biāo)，但必須再經(jīng)過計算才能獲得二維參數(shù)空間的值。因此，該類單元不能從現(xiàn)有的有限元法前處理的商業(yè)軟件中直接得到。為了保證良好的積分和插值精度，一般要求二維參數(shù)單元在對應(yīng)的三維物理空間中的網(wǎng)格具有良好的形狀。為此，作者將鋪磚法運用到任意復(fù)雜曲面生成參數(shù)曲面網(wǎng)格，實現(xiàn)邊界面法的前處理。

UG是一款采用全參數(shù)化建模CAD軟件，并提供了強大的二次開發(fā)工具UG/Open[8]。本文通過VC++編程調(diào)用UG/Open庫函數(shù)獲取CAD模型的曲面參數(shù)信息。在參數(shù)域內(nèi)利用改進(jìn)的鋪磚法生成網(wǎng)格，然后將二維參數(shù)域上生成的網(wǎng)格映射到曲面的三維物理空間，從而生成曲面的最終網(wǎng)格。為了使二維參數(shù)域上生成的網(wǎng)格投到曲面的三維物理空間保持較好的質(zhì)量，本文在參數(shù)域上生成網(wǎng)格時運用黎曼度量方法對網(wǎng)格的形狀和大小進(jìn)行控制。文章最后給出了4個算例驗證了算法的可靠性，并且說明了不論是在平面還是曲面上都能生成質(zhì)量較好的四邊形網(wǎng)格。

1 黎曼度量

曲面二維參數(shù)域上的點與三維物理空間的點是一一對應(yīng)的關(guān)系。對于簡單的曲面直接映射法能生成質(zhì)量較高的網(wǎng)格，對于曲率變化劇烈的曲面，由于其映射函數(shù)在參數(shù)域的兩個方向的非線性生成的網(wǎng)格的質(zhì)量往往較差。為了使二維參數(shù)域上生成的網(wǎng)格投影到曲面的三維物理空間保持較好的質(zhì)量，采用黎曼度量的方法可以控制網(wǎng)格形狀和大小。

1.1 參數(shù)曲面的黎曼度量

三維空間的一個參數(shù)曲面S與二維參數(shù)空間的一一對應(yīng)關(guān)系可用函數(shù)表示為

在UG二次開發(fā)環(huán)境中，只需知道二維參數(shù)域中每個節(jié)點u和v兩個參數(shù)值就可以調(diào)用的庫函數(shù)UF_MODL_evaluate_face（…）獲得節(jié)點沿參數(shù)u和v的切向量Su和Sv，進(jìn)而得到參數(shù)曲面黎曼度量矩陣。

1.2 參數(shù)線段弧長計算

參數(shù)域內(nèi)線段兩端點為 P1和 P2，則該線段上任意點可表示為 P = P1+ ( P2- P1)? t ，0≤t≤1。P點的黎曼度量可表示為

考慮黎曼度量后線段上兩點A和B的在三維物理空間長度可以用下面公式計算[2]

當(dāng)曲線上兩點A( t1)，B( t2)的度量矩陣特征值變化不大時可用上式，否則需要把線段分成若干段，計算每小段li的長度，總長度為

1.3 參數(shù)向量的旋轉(zhuǎn)

二維參數(shù)空間中向量AB需逆時針旋轉(zhuǎn)一個角度后得到新的向量AC，則向量AC由下式計算得到[3-4]

lAB為向量AB在三維物理空間的長度，可根據(jù)參數(shù)弧長計算公式(3a)或者式(3c)得到，lideal為生成C點所需的三維物理空間理想長度，φ根據(jù)前沿DA和AB映射到三維物理空間中形成的角度∠DAB而確定。如圖1(a)圖1(b)，沿AB逆時針方AB逆時針方向φ為的度量矩陣，根據(jù)公式(4)即可求出向量AC，點A為已知點，則可以求出C點在二維參數(shù)空間的值。

圖1 向量旋轉(zhuǎn)角度

2 曲面參數(shù)空間四邊形生成

2.1 網(wǎng)格生成過程

根據(jù)鋪磚法的一般原理，如圖2給出了生成曲面網(wǎng)格的整個流程。本文對相交處理做了改進(jìn)，并且在2.2做了詳細(xì)的闡述，其它過程只簡單介紹，其細(xì)節(jié)可以參考文獻(xiàn)[5, 9-10]。

1）邊界離散

調(diào)用UG庫函數(shù)獲得實體表面的參數(shù)域，對參數(shù)域的邊界進(jìn)行離散。內(nèi)邊界點按順時針方向連接，外邊界點按逆時針方向，且保證離散點數(shù)為偶數(shù)。

2）節(jié)點分類及最優(yōu)邊選擇

節(jié)點內(nèi)角是指節(jié)點Ni與其所在邊界上前一節(jié)點Ni-1和后一節(jié)點Ni+1投到三維物理空間的角度。根據(jù)節(jié)點內(nèi)角的不同分為 4種類型：(1) 終 止 節(jié) 點 α≤ 1 35°；(2) 邊 節(jié) 點135°＜ α ≤ 2 25°；(3) 角節(jié)點225° ＜ α ≤ 3 00°；(4)轉(zhuǎn)節(jié)點300° ＜ α ≤ 3 60°。

如果Ni為終止節(jié)點則邊界 NiNi+1為最優(yōu)邊開始生成單元，每次選取前沿邊界中的一條邊作為基邊生成四邊形。

3）生成新單元及相交判斷

本文中單元是在二維參數(shù)域中生成的。每次生成一個新的單元是根據(jù)基邊兩個端點的節(jié)點內(nèi)角在參數(shù)域中生成新的節(jié)點，從而生成新的單元。新節(jié)點的位置可根據(jù)公式(4)確定，其中 lAB為基邊在三維域中長度，lideal的長度可以參考文獻(xiàn)[5, 9-10]中計算新節(jié)點的向量長度，MA是節(jié)點的度量矩陣。形成單元后，新邊要判斷是否與前沿邊界相交，如果相交就進(jìn)行處理，在下節(jié)中詳細(xì)介紹相交處理過程。

4）生成一排單元及邊界調(diào)整

以終止節(jié)點作為新生單元的起點，如果沒有終止節(jié)點選擇最小節(jié)點內(nèi)角作為起點，直到再次遇到終止節(jié)點完成一排單元的生成，當(dāng)所有的活動前沿邊界都作為基邊生成單元后就完成一層單元的生成，然后更新活動前沿邊界。一層單元生成后，如果更新的活動前沿邊界尺寸有變大或者縮小的趨勢則進(jìn)行單元楔入或者刪除單元操作，具體過程參考文獻(xiàn)[5, 9-10]。

5）光滑處理和縫合處理

為了使生成的邊界尺寸大小均勻，以及邊的垂直度，可以對更新的前沿邊界節(jié)點進(jìn)行邊界光滑，內(nèi)部節(jié)點進(jìn)行內(nèi)部光滑。如果尺寸過渡比較大或者出現(xiàn)內(nèi)角很小的情況可以進(jìn)行縫合處理，詳細(xì)情況可以參考文獻(xiàn)[5, 9-10]。

6）閉合處理及優(yōu)化網(wǎng)格

當(dāng)前沿邊界節(jié)點數(shù)目小于等于4時就進(jìn)行封閉處理，整個區(qū)域網(wǎng)格生成完成。用改進(jìn)的鋪路法完成全四邊形網(wǎng)格生成后，每個節(jié)點最理想的是由4個單元共用，這樣形成的四邊形較規(guī)則。如果每個節(jié)點由3～5個單元共用，生成的單元是符合要求的。由于幾何邊界的復(fù)雜性和網(wǎng)格尺寸不均勻常常導(dǎo)致局部網(wǎng)格形狀不好，所以需進(jìn)一步優(yōu)化。優(yōu)化后需再對除幾何邊界節(jié)點外節(jié)點用內(nèi)部節(jié)點光滑方法進(jìn)行幾次迭代光滑，保證生成的網(wǎng)格質(zhì)量更好。

圖2 程序流程

2.2 相交處理

原有的鋪磚法網(wǎng)格是逐排生成的，生成一排單元后再去處理相交問題。如果采用這種方法處理，特別是出現(xiàn)多次相交還要進(jìn)行優(yōu)先處理判斷，這樣使程序比較復(fù)雜。本文對相交處理過程做了改進(jìn)，采用生成一個單元時就進(jìn)行相交判斷，如果遇到相交就進(jìn)行相交處理。

經(jīng)過改進(jìn)后的相交處理過程可分為 3個部分：

1）優(yōu)先處理新單元中相交線段新生成的點；

2）如果1）未完成處理則處理新單元中相交線段原有節(jié)點；

3）如果前兩步都未處理好，則搜索基邊節(jié)點附近的點，如果滿足要求則替代新節(jié)點，如果沒有滿足要求的點則取消此單元生成，并以相鄰邊作為基邊生成單元。

如圖 3所示，新生成的單元邊 NiNi+1和NiNi-1與前沿ab相交，Ni為新生成的點，相交處理后必須保證邊界形成環(huán)后的節(jié)點數(shù)為偶數(shù)，圖3中是兩個環(huán)合并為一個環(huán)，Ni點距離b點近所以選b作為替代的點，圖3(c)所示。點Ni被點b替代后根據(jù)點b兩邊的角度來確定邊的合并，如圖3(d)所示為合并邊后的結(jié)果。

圖3 新節(jié)點相交處理

如圖4所示，圖4(a)中直線未相交，此時以邊NiNi-1作為基邊生成四邊形，生成新節(jié)點Nj并生成四邊形，圖4(b)圖所示，新生成的單元邊Ni+1Nj與前沿ab相交，如果新節(jié)點由節(jié)點a替代形成的邊界環(huán)不滿足偶數(shù)個節(jié)點要求，如果新節(jié)點由節(jié)點b替代生成的邊還是與前沿相交，此時考慮節(jié)點 Ni+1由a或者b替代。在本例中由a替代，如圖4(c)所示?？p合最終生成的圖形如圖4(d)所示。

圖4 原節(jié)點相交處理

如果以上兩個處理過程都不能生成單元，則進(jìn)行第3步處理。如圖5所示，邊界AB為基邊，根據(jù)節(jié)點B的節(jié)點內(nèi)角生成C點，假設(shè)C點按前面兩步都不能處理相交問題，則以B為中心在參數(shù)空間建立一個搜索區(qū)域，在區(qū)域內(nèi)尋找最優(yōu)的點。

首先確定在三維空間的搜索半徑r3D，一般取邊界AB在三維空間長度的兩倍?？紤]度量矩陣，以搜索半徑r3D構(gòu)成的區(qū)域在參數(shù)空間里為一橢圓，參數(shù)橢圓的中心在B點。因此，要計算參數(shù)橢圓的長半徑和短半徑，以確定簡單的參數(shù)矩形。根據(jù)搜索半徑r3D，由以下公式計算對應(yīng)的參數(shù)半徑

方程(5)等價于中心在圓點極坐標(biāo)中的橢圓方程，能轉(zhuǎn)化為參數(shù)坐標(biāo)下橢圓形式的方程為

通過旋轉(zhuǎn)橢圓與坐標(biāo)軸平行，方程(6)可以簡化，得到橢圓在參數(shù)空間的長短半徑分別定義為

得到橢圓的長短半徑后，可以簡單地在參數(shù)空間里確定矩形搜索區(qū)域，如圖4所示的矩形框。搜索區(qū)域建立好后選擇最優(yōu)點，如圖6所示，前沿節(jié)點P1, P2, P3為搜索區(qū)域內(nèi)的點，最優(yōu)點形成四邊形后要保證每個邊界前沿節(jié)點數(shù)目為偶數(shù)個，圖中P2不滿足要求，P1和P3滿足奇偶性要求，本文中選取∠APB張角最大的點，即P1是最優(yōu)點。

圖6 選擇最優(yōu)點

3 四邊形質(zhì)量評價因子

網(wǎng)格質(zhì)量的好壞將影響計算精度，本文采用幾何準(zhǔn)則對四邊形網(wǎng)格進(jìn)行質(zhì)量評價。如圖7所示四邊形ABCD的兩條對角線將其分為4個三角形△ABC，△ABD，△CDB，△CDA，通過計算這4個三角形的質(zhì)量因子得出該四邊形的質(zhì)量因子。

圖7 分割四邊形

三角形ABC的質(zhì)量α因子[11]定義如下

n是三角形的單位法向量，α的意義就是三角形面積與變長的平方和之比，0＜α≤1，α越大三角形質(zhì)量越好。令4個三角形△ABC，△ABD，△CDB，△CDA的α因子分別為1α，2α，3α，四邊形β因子[12]定義為

0＜β≤1，β越大四邊形質(zhì)量越好，根據(jù)文獻(xiàn)[11]中β不同的范圍來評價四邊形質(zhì)量，如表1所示。

表1 四邊形 β 評價因子[11]

4 UG實現(xiàn)過程

本文利用VC2005編程平臺通過External連接方式實現(xiàn)UG二次開發(fā)。首先，在VC2005環(huán)境里建立一個動態(tài)(.dll)項目。然后，對項目進(jìn)行設(shè)置即可與 UG/Open API建立連接關(guān)系。在VC2005編譯環(huán)境里對VC++項目的設(shè)置為：

1）菜單選項-工具選項項目和解決方案VC++目錄：

在“庫文件”和“包含文件”里添加UG安裝目錄下的UGSNX 4.0UGOPEN

2）菜單選項-項目屬性頁連接器輸入附加依賴項：

添 UG Open 的庫文件名 libufun. lib,libugopenint. lib和libnxopencpp.lib。

3）菜單選項-項目屬性頁連接器常規(guī)輸出文件：

添加‥startup$(ProjectName).dll。

4）菜單選項-項目屬性頁配置屬性常規(guī)字符集：

下拉選項中選擇為“使用多字節(jié)字符集”。

通過上面的設(shè)置，VC++環(huán)境里的項目與UG/Open建立了聯(lián)系，成功編譯項目后就在該項目目錄下的startup文件中生成相應(yīng)可以在UG中執(zhí)行的.dll文件。在 UG環(huán)境里利用菜單工具——調(diào)用.dll就可以在相應(yīng)的曲面上生成網(wǎng)格。

當(dāng)與 UG/Open建立聯(lián)系后可以調(diào)用其中的庫函數(shù)獲取實體模型的曲面信息，根據(jù)獲得的曲面參數(shù)空間進(jìn)行邊界離散，并運用第2節(jié)中的網(wǎng)格生成方法在參數(shù)空間生成網(wǎng)格，并映射到三維物理空間。

下面給出在 UG/Open API中常訪問和生成幾何的函數(shù)以及功能：

UF_MODL_ask_body_faces(…)：獲得實體的所有曲面；

UF_MODL_ask_face_edges(…)：獲得曲面的所有邊界曲線；

UF_MODL_ask_face_uv_minmax(…)：獲得曲面的參數(shù)范圍；

UF_MODL_evaluate_face（…）：獲得節(jié)點沿參數(shù)u和v的切向量以及節(jié)點映射到三維物理空間的坐標(biāo)；

UF_POINT_create_on_surface(…)：在曲面上生成點；

UF_CURVE_create_line(…)：生成線段；

UF_OBJ_set_color(…)：生成的網(wǎng)格在物理空間的顏色顯示。

5 網(wǎng)格生成實例

該節(jié)給出了利用本文算法生成曲面四邊形網(wǎng)格實例，如圖8～圖11所示。網(wǎng)格單元直接在UG環(huán)境中的CAD模型表面生成，保持了模型的原始幾何信息。根據(jù)網(wǎng)格質(zhì)量因子對各種實例的網(wǎng)格質(zhì)量進(jìn)行了評價，其結(jié)果如表2所示?？梢钥闯龌诶杪攘扛倪M(jìn)后的鋪磚法無論在平面還是曲率大的曲面都可以獲得質(zhì)量較好的網(wǎng)格。

圖8 實例1

圖9 實例2

圖10 實例3

圖11 實例4

表2 實例質(zhì)量評價表

6 結(jié)束語

本文基于黎曼度量將鋪磚法擴展運用到任意曲面的四邊形網(wǎng)格生成，并在原有算法的基礎(chǔ)上對相交處理過程做了改進(jìn)，使得相交處理更加簡單有效?；赨G二次開發(fā)工具，利用VC2005平臺在UG建模環(huán)境中實現(xiàn)了曲面的二維參數(shù)空間網(wǎng)格劃分，且網(wǎng)格單元在三維空間質(zhì)量較好。在網(wǎng)格生成過程中，調(diào)用UG Open/API中相應(yīng)的庫函數(shù)，獲取CAD實體邊界曲面參數(shù)信息，并將最終生成的網(wǎng)格顯示在實體表面。后續(xù)的工作是將這種參數(shù)四邊形網(wǎng)格應(yīng)用到邊界面法中作為前處理網(wǎng)格。邊界面法與CAD無縫連接能避免幾何誤差，提高計算精度，實現(xiàn)CAD與CAE的無縫連接，對數(shù)值計算有非常重要的意義。

[1]張建華, 葉尚輝. 有限元網(wǎng)格自動生成典型方法及發(fā)展方向[J]. 計算機輔助設(shè)計與制造, 1996, (2):28-31.

[2]關(guān)振群, 單菊林, 顧元憲. 基于黎曼度量的復(fù)雜參數(shù)曲面有限元網(wǎng)格生成方法[J]. 計算機學(xué)報, 2006,29(10): 1823-1833.

[3]Lee C K. Automatic mesh generation using metric advancing front approach [J]. Engineering Computations, 1999, (16): 230-263.

[4]Lee C K. Automatic adaptive metric advancing front triangulation over curved surfaces [J]. Engineering Computations, 2000, 17: 48-74.

[5]Blacker T D, Stephenson M B. Paving: a new approach to automated quadrilateral mesh generation [J].Numer. Meth. Engrg, 1991, 32: 811-847.

[6]Zhang Jianming, Qin Xianyun, Xu Han, et al. A boundary face method for potential problems in three dimensions [J]. Numer. Meth. Engrg, 2009, 80:320-337.

[7]覃先云. 基于參數(shù)曲面的邊界面法研究[D]. 長沙:湖南大學(xué), 2009.

[8]侯永濤, 丁向陽. UG/Open 二次開發(fā)與實例精解[M].北京: 化學(xué)工業(yè)出版社, 2002: 1-2.

[9]張清萍, 尚 勇, 趙國群. 二維全四邊形網(wǎng)格的自動生成算法[J]. 山東大學(xué)學(xué)報, 2002, 32(3):254-259.

[10]林勝良, 方 興, 張 武, 等. 一種改進(jìn)的高品質(zhì)全四邊形網(wǎng)格生成方法[J]. 江南大學(xué)學(xué)報, 2006,5(1): 70-73.

[11]Lee C K, Lo S H. A new scheme for the generation of a graded quadrilateral mesh [J]. Computers and Structures, 1994, 52: 847-857.