汪 攀 張見明 韓 磊 鞠傳明

池寶濤湖南大學汽車車身先進設計制造國家重點實驗室，長沙，410082

0 引言

網(wǎng)格生成是數(shù)值模擬的主要性能瓶頸，其自動生成算法的研究一直被廣泛關注。有限元方法的成功應用在很大程度上取決于對分析對象進行有限元網(wǎng)格劃分的合理性。相對于四面體網(wǎng)格，六面體網(wǎng)格在計算精度、劃分數(shù)量、抗畸變程度以及重劃分次數(shù)等方面均具有優(yōu)勢[1]，因此六面體網(wǎng)格也稱為黃金網(wǎng)格。但是，由于六面體網(wǎng)格復雜的拓撲關系以及較差的幾何自適應能力，故針對復雜三維實體的全六面體單元網(wǎng)格全自動生成方法，目前仍處于探索階段。目前有代表性的全六面體網(wǎng)格自動生成方法有：超限映射法[2]、掃描法[3-4]、基于柵格法[5-7]、前沿推進法[8-9]和多子區(qū)域法[10]。其中超限映射法的優(yōu)點是算法簡單、速度快、單元質(zhì)量好、密度可控制，并且可與形狀優(yōu)化算法集成，因此，映射法在眾多的商業(yè)有限元分析軟件中占有重要的地位。但是，映射法一般只能直接處理簡單的單連通域問題，并且對于帶有倒圓角、小孔等小特征的三維實體，其生成的六面體網(wǎng)格質(zhì)量不好；對于復雜的多連通域問題，通常需要首先手工將待剖分域分解成幾何形狀規(guī)則的可映射子區(qū)域，然后在每個子區(qū)域上應用超限映射法，這一缺陷不利于大規(guī)模全自動網(wǎng)格劃分的實現(xiàn)。為了解決這一問題，PRICE等[11]提出用中面法將三維待剖分域分解成簡單子區(qū)域，但是現(xiàn)有中面算法一般需要進行大量幾何與代數(shù)計算，自動化程度和幾何適應能力也有待于提高，且對于存在多凹面的幾何體來說中軸面生成存在困難。

針對上述問題，本文提出了一種基于子域分解的全六面體網(wǎng)格自動生成方法，該方法首先提取實體模型的分解特征，然后通過分解特征形成分解面，最后利用分解面將復雜實體分解為多個可映射的簡單子域，并在各子域上用超限映射法生成六面體網(wǎng)格。

1 算法原理

1.1 相關定義

1.1.1 與邊相關的定義

(1)二面角：與邊相連的兩個面的夾角,設為α。

(2)凸邊：當π/2-ε<α<π/2+ε時，則該邊為凸邊，其中ε取π/6，下同。

(3)自然邊：當π-ε<α<π+ε時，則該邊為自然邊。

(4)凹邊：當3π/2-ε<α<3π/2+ε時，則該邊為凹邊。

(5)反轉(zhuǎn)邊：當3π/2+ε<α時，則該邊為反轉(zhuǎn)邊。

(6)特征邊：凹邊和反轉(zhuǎn)邊統(tǒng)稱為特征邊。

(7)虛邊：用于形成封閉的分解環(huán)時構(gòu)成出來的邊。

1.1.2 與環(huán)相關的定義

(1)特征環(huán)：由相連的特征邊構(gòu)成的環(huán)。特征環(huán)有可能是開環(huán)，也有可能是封閉的環(huán)。

(2)分解環(huán)：由特征環(huán)和虛邊組成的連續(xù)封閉的環(huán)。用于形成分割面，對三維實體進行分割。

1.1.3 與面相關的定義

分解面：由分解環(huán)得到的虛擬面，用于分割三維實體。

1.2 算法流程

本文所提出算法的主要步驟如下。

(1)識別實體的小特征，如圓倒角、同心圓、圓孔等，得到其等效的幾何模型。

(2)從三維實體的等效模型中提取實體的特征邊(凹邊和反轉(zhuǎn)邊)，并將其加入特征邊鏈表。

(3)遍歷特征邊鏈表，找到相互連接的特征邊，構(gòu)成一個特征環(huán)，并將其加入特征環(huán)鏈表。注意：特征環(huán)有可能包含多連邊，也有可能只有一條邊，有可能是連續(xù)封閉的環(huán)，也有可能不是封閉的環(huán)。

(4)遍歷特征環(huán)鏈表。若取到的特征環(huán)是連續(xù)封閉的環(huán)，則直接通過該特征環(huán)及其支撐面構(gòu)造出分解面；若取到的特征環(huán)不是封閉的，則需要構(gòu)造合適的虛邊，將特征環(huán)補成連續(xù)封閉的環(huán)，然后再構(gòu)造分解面。

(5)通過分解面將三維實體分割為兩個部分，為了保證最終生成的體網(wǎng)格連續(xù)一致，需要建立這兩部分實體與分解面之間的拓撲關系。

(6)對分割后的實體進行網(wǎng)格劃分時，為了保證空間網(wǎng)格的一致性，要先在分解面上生成面網(wǎng)格，再映射到被其分割的實體表面，之后用超限映射生成六面體網(wǎng)格。

(7)最后將各子體的六面體網(wǎng)格數(shù)據(jù)合并，即得到整個三維實體的網(wǎng)格數(shù)據(jù)。

2 實現(xiàn)過程

2.1 處理小特征

2.1.1 識別小特征

對于很多機械產(chǎn)品，為了避免應力集中，都會引入倒圓角設計，這會導致捕獲實體模型的特征邊失敗。如圖1a所示的邊a、b，其二面角α=π，此時會將a、b兩條邊鑒別為自然邊，但是，很明顯圖1a所示的模型和圖1b是等效的，此時體分解失敗，所以需對倒圓角作特殊處理。要處理倒圓角，首先要識別倒圓角面，下面給出倒圓角面的三個條件：①面包含四條邊；②有兩條圓弧邊、兩條直線邊；③兩條圓弧邊是對邊，且其對應的圓心角相同。若給定的實體表面滿足以上三個條件，則可以判定為倒圓角面。

(a)原始模型 (b)等效模型 圖1 倒圓角的處理Fig.1 The treatment of fillet

2.1.2 模板

機械產(chǎn)品中，經(jīng)常會遇見帶圓倒角、同心圓柱面、圓孔等小特征的實體，這些小特征的存在阻礙了六面體網(wǎng)格的生成，或使得生成的六面體網(wǎng)格質(zhì)量不好，為了解決此問題，本文給出了如圖2～圖4所示的模板，圖中所標識的點在三維空間是一條直線，引入模板以后，可以準確地識別邊的類型，從而可以有效合理地進行體分解流程。其中，圖2所示的模板可以準確地提取三維實體模型的特征邊，從而對三維實體進行有效體分隔，圖3所示的同心圓模板和圖4所示的圓孔模板，可以顯著地提高生成六面體網(wǎng)格的質(zhì)量。

圖2 倒圓角的模板Fig.2 The templates of fillet

圖3 同心圓的模板Fig.3 The templates of concentric circles

(a)幾何模型 (b)分解模型 (c)網(wǎng)格劃分圖4 圓孔的模板Fig.4 The templates of hole

2.2 形成合理的分解面

2.2.1 提取特征邊

特征邊是三維實體上形成分解環(huán)的邊，本文方法中所指的分解邊是指凹邊和反轉(zhuǎn)邊。凹邊和反轉(zhuǎn)邊是形成分解環(huán)，最終形成分解面的基礎，準確獲取實體模型的分解邊是進行體分解的關鍵步驟。獲取三維實體的分解邊主要包括以下兩個步驟：

(1)識別三維實體的小特征，根據(jù)圖2～圖4所示的模板，得到其等效模型。

(2)遍歷等效模型的邊界邊，計算其對應的二面角，并標記其類型，將凹邊和反轉(zhuǎn)邊存在鏈表中，如圖5所示的粗線邊。

圖5 圓倒角的特征邊提取示意圖Fig.5 The feature abstraction of fillet

2.2.2 形成分解環(huán)

獲取了三維實體的特征邊以后，遍歷得到的特征邊，相互連接的特征邊組成一個分解環(huán)。當分解環(huán)是閉環(huán)時，可直接由分解環(huán)形成分解面，然后對三維實體進行有效分割；當分解環(huán)是開環(huán)時，需通過搜索形成合適的“虛擬邊”，將開環(huán)補成連續(xù)封閉的分解環(huán)。封閉的分解環(huán)才能有效地對三維實體進行分割。如圖6a所示的實體模型，根據(jù)上述邊的定義提取到四條特征邊(凹邊)，如圖6b所示的粗線段。這4條特征邊首尾相連，構(gòu)成一個封閉的分解環(huán)，可直接形成分解面對三維實體進行分割。

(a) (b) 圖6 分解環(huán)形成示意圖(一)Fig.6 The generation of parting loop

從圖5所示模型的等效模型中獲取的特征邊如圖7a所示的粗線段AB，通過該分解邊形成特征環(huán)時，沒有搜索到與其相連的分解邊，此時該特征邊形成的特征環(huán)即為一條孤立的線段，此時需要構(gòu)造虛邊。具體方法如下：在特征邊的端點處，作與其相連的兩條邊的延長線，并與面上另一條邊相交，該端點與交點即為一條虛邊，如7b所示的粗線段AC和AD，重復上述過程，直到構(gòu)建的虛邊可以和分解邊構(gòu)成一個封閉的環(huán)。特征邊沿著兩個方向構(gòu)造虛邊，可以得到兩個分解環(huán)，如圖7c所示的ADEB和圖7d所示的ACFB，此時需要選擇一個相對合適的分解環(huán)。選取的準則如下：設分解環(huán)上最短邊長度為Lmin，最長邊長度為Lmax，比例β=Lmin/Lmax，取β較大的分解環(huán)，β越大，分解環(huán)的形狀越趨向正方形(當分解環(huán)是四條邊的時候)，合理地選擇分解環(huán)可以有效地避免分割后的實體出現(xiàn)狹長的片體，從而不利于生成高質(zhì)量的六面體網(wǎng)格。

(a) (b)

(c) (d)圖7 分解環(huán)形成示意圖(二)Fig.7 The generation of parting loop

2.2.3 形成分解面

三維實體之間是通過面連接的，所以需要由面來對實體進行分割。形成分解環(huán)以后，如果分解環(huán)中存在虛邊，則根據(jù)分解環(huán)構(gòu)建一個虛平面；如果分割環(huán)中不存在虛邊，則先獲取分割環(huán)在三維實體上的支撐表面(該支撐表面不一定是平面)，得到支撐表面以后，根據(jù)支撐表面和分解環(huán)構(gòu)造出虛面，此時的構(gòu)造并不是真的生成一個面，而是在支撐面的基礎上加了一個約束邊界，這個虛面的幾何是基于支撐面，拓撲是基于分解環(huán)。分解面形成如圖8所示。

圖8 分解面形成示意圖Fig.8 The generation of parting plane

2.3 體分解

對于復雜實體，可能會形成多個分解面，如何利用這些分解面對三維實體進行全自動分割是本文算法的關鍵。為了簡化程序?qū)崿F(xiàn)過程，本文將所有分解面存儲在鏈表中，然后遍歷該鏈表。對于某一分解面，將目標域分解為兩個子域，若這兩個子域中，至少有一個可以利用超限映射法進行六面體網(wǎng)格劃分(子域包含6個實體面或為圖3、圖4所示的模板)，則該分解過程合法，即可以利用該分解面對目標域進行分解，否則將該分解過程延后至下一次循環(huán)。上述操作可以形象地描述為，從復雜的目標域上，逐步分割出一個可用映射法進行網(wǎng)格劃分的子域。利用上述方法進行有序的分解，可以避免遞歸算法，從而節(jié)省了內(nèi)存開銷，且方便程序的調(diào)試與維護。分解面將目標域分解為兩個子域以后，需要建立兩個子域與分解面的拓撲關系。目前，大多數(shù)CAD軟件采用基于邊界表征[12]的數(shù)據(jù)結(jié)構(gòu)(B-rep)。在邊界表征的數(shù)據(jù)結(jié)構(gòu)中，實體(body)是由面(face)表示的，面(face)是由環(huán)(loop)表示的，環(huán)(loop)是由邊(edge)表示的。本文依據(jù)邊界表征的數(shù)據(jù)結(jié)構(gòu)構(gòu)造分解面與兩個子域的拓撲關系，首先由分解邊構(gòu)造分解面，然后將分解面存儲于子域中用于存儲幾何面數(shù)組FaceArray。

2.4 子域網(wǎng)格劃分

2.4.1 子域邊界離散

用超限映射法生成六面體時，要求三維實體的幾何邊界離散線段數(shù)相等，因此對于每個分割后的塊體，在每個面上其對應的離散線段數(shù)要相等。這實際上就是一個帶約束的線性規(guī)劃問題，通過求解下列方程，即可得到所有邊界的離散線段數(shù)。

目標函數(shù)為

約束條件為

式中，k為進行體分解以后，三維實體包含的邊；Li為第i條邊的實際長度；Ls為自適應的網(wǎng)格尺寸值；iFX+為面上方向與X軸正向保持一致的所有邊界的集合；iFX-為面上方向與X軸負向保持一致的所有邊界的集合；iFY+為方向與Y軸正向保持一致的所有邊界的集合；iFY-為方向與Y軸負向保持一致的所有邊界的集合。

上述問題可以采用整數(shù)規(guī)劃方法來求解[13]。

2.4.2 超限映射法

對三維實體進行分割以后得到的可映射子體，可以用超限映射法[13]生成六面體網(wǎng)格。對于任意的包含6個實體表面，且是單連通域的三維實體，其對應的參數(shù)空間為[0,1]×[0,1]×[0,1]，設(ξ,η,γ)為參數(shù)空間[0,1]×[0,1]×[0,1]中任意的一點，則其對應在物理空間的坐標是由(ξ,η,γ)在6個表面上對應的坐標和8個頂點決定的[14]。對于任意的包含4條邊界，且是單連通域的面，其對應的參數(shù)空間為[0,1]×[0,1]，設(u,v)為參數(shù)空間[0,1]×[0,1]中任意的一點，則其對應在物理空間的坐標是由(u,v)在4條邊界上對應的坐標和“四邊形”面的4個頂點決定的[13]，現(xiàn)以面的超限映射法加以說明。設“四邊形”面的參數(shù)域如圖9a所示，則(u,v)對應的物理域坐標可由下式得到：

式中，f1、f2、f3、f4為實體表面上四條邊界，如圖9b所示。

(a)參數(shù)域

(b)物理域圖9 “四邊形”面的物理域及參數(shù)域Fig.9 The physical domain and parameter domain of four sides face

2.5 子域網(wǎng)格合并

在各子域上使用超限映射法進行網(wǎng)格劃分以后，需要將子域網(wǎng)格合并，從而得到整個目標域的網(wǎng)格。本文首先對分解面進行網(wǎng)格劃分，并將其作為包含該分解面的兩個子域在該表面上的網(wǎng)格數(shù)據(jù);其次，對子域上其他表面進行網(wǎng)格劃分，并將生成的表面網(wǎng)格數(shù)據(jù)加入存儲目標域網(wǎng)格數(shù)據(jù)的容器中；然后，依次在各子域上采用超限映射法生成域內(nèi)的網(wǎng)格數(shù)據(jù)，并將生成的域內(nèi)網(wǎng)格數(shù)據(jù)加入目標域網(wǎng)格數(shù)據(jù)容器中。按照上述方式得到的網(wǎng)格數(shù)據(jù)在分解面上是連續(xù)一致的，且界面處的網(wǎng)格數(shù)據(jù)只存儲一遍。界面上的網(wǎng)格生成示意圖見圖10。

圖10 界面上的面網(wǎng)格生成Fig.10 The mesh generation of interface

3 數(shù)值算例

下面給出三個機械零部件的六面體網(wǎng)格劃分實例，如圖11～圖13所示。三個實體模型不能直接利用商業(yè)軟件對其進行六面體劃分，必須要進行人工分割，分解為簡單子域才能有效生成六面體網(wǎng)格。本文利用特征識別技術可以全自動地進行體分解，得到一系列簡單規(guī)則可映射的子域，然后利用超限映射法生成質(zhì)量較好的結(jié)構(gòu)化六面體網(wǎng)格，實現(xiàn)自動劃分。

(a)原始模型

(b)分解模型

(c)生成的六面體網(wǎng)格圖11 六面體網(wǎng)格生成算例一Fig.11 The first example of hexahedral mesh generation

(a)原始模型 (b)分解體

(c)進一步分解體 (d)生成的六面體網(wǎng)格圖12 六面體網(wǎng)格生成算例二Fig.12 The second example of hexahedral mesh generation

圖13 六面體網(wǎng)格生成算例三Fig.13 The third example of hexahedral mesh generation

4 總結(jié)

(1)由于子域網(wǎng)格是由超限映射法生成的，因此本文提出的方法能夠生成質(zhì)量較好的結(jié)構(gòu)化六面體網(wǎng)格，且算法效率較高。

(2)對于含有倒角、同心圓、小孔等小特征的三維實體，本文引入的模板能夠生成質(zhì)量較好的六面體網(wǎng)格。

(3)相對于傳統(tǒng)的商業(yè)軟件，本文的算法能夠全自動地生成六面體網(wǎng)格，避免了大量人機交互操作，提高了網(wǎng)格生成的效率。

(4)數(shù)值實驗表明，對于無法用商業(yè)軟件直接進六面體劃分的復雜實體，本文提出的方法能夠全自動地進行六面體網(wǎng)格劃分。

[1] SCHONNING A, OOMMEN B, IONESCU I, et al. Hexahedral Mesh Development of Free-formed Geometry: the Human Femur Exemplified[J]. Computer-aided Design,2009,41(8):566-572.

[2] GORDON W J, THIEL L C. Transfinite Mappings and Their Application to Grid Generation[J]. Applied Mathematics & Computation,1982,10/11(4):171-233.

[3] MUKHERJEE N, PEDDI B, CABELLO J, et al. Automatic Hexahedral Sweep Mesh Generation of Open Volumes[C]//International Meshing Roundtable. Orlando,2013:333-347.

[4] 陳建軍, 肖周芳, 曹建,等. 多源掃掠體全六面體網(wǎng)格自動生成算法[J].浙江大學學報工學版,2012,46(2):274-279. CHEN Jianjun, XIAO Fangzhou, CAO Jian, et al. Automatic Hexahedral Mesh Generation for Many-to-one Sweep Volume[J]. Journal of Zhejiang University,2012,46(2):274-279.

[5] ZHANG H, ZHAO G. Adaptive Hexahedral Mesh Generation Based on Local Domain Curvature and Thickness Using a Modified Grid-based Method[J]. Finite Elements in Analysis & Design,2007,43(9):691-704.

[6] 黃麗麗,趙國群. 基于柵格法的三維六面體網(wǎng)格質(zhì)量優(yōu)化[J]. 中國機械工程,2009,20(21):2603-2608. HUANG Lili, ZHAO Guoqun. Optimization of Grid-based Three-dimensional Hexahedral Meshes[J]. China Mechanical Engineering,2009,20(21):2603-2608.

[7] SUN L, ZHAO G. Adaptive Hexahedral Mesh Generation and Quality Optimization for Solid Models with Thin Features Using a Grid-based Method[J]. Engineering with Computers,2016,32(1):61-84.

[8] KREMER M, BOMMES D, LIM I, et al. Advanced Automatic Hexahedral Mesh Generation from Surface Quad Meshes[C]// Proceedings of International Meshing Roundtable.Orlando,2013:147-164.

[9] KAWAMURA Y, ISLAM M S, SUMI Y. A Strategy of Automatic Hexahedral Mesh Generation by Using an Improved Whisker-weaving Method with a Surface Mesh Modification Procedure[J]. Engineering with Computers,2008,24(3):215-229.

[10] TAM T K H, ARMSTRONG C G. Finite Element Meshcontrol by Integer Programming[J]. International Journal for Numerical Methods in Engineering,1993,36(15):2581-2605.

[11] PRICE M A, ARMSTRONG C G, SABIN M A. Hexahedral Mesh Generation by Medial Surface Subdivision: Part I. Solids with Convex Edges[J]. International Journal for Numerical Methods in Engineering,1995,38(19):3335-3359.

[12] 孫家廣. 計算機圖形學[M]. 3版.北京: 清華大學出版社, 1998. SUN Jiaguang. Computer Graphics[M]. 3rd ed. Beijing: Tsinghua University Press,1998.

[13] WHITELEY M, WHITE D, BENZLEY S, et al. Two and Three-quarter Dimensional Meshing Facilitators[J]. Engineering with Computers,1996,12(3):144-154.

[14] LI T S, ARMSTRONG C G, MCKEAG R M. Quad Mesh Generation fork-sided Faces and Hex Mesh Generation for Trivalent Polyhedral[J]. Finite Elements in Analysis & Design,1997,26(4):279-301.