李海峰 吳冀川 劉建波 梁宇兵

1.中國工程物理研究院計算機應(yīng)用研究所，綿陽，6219002.新加坡國立大學，新加坡，119077

0 引言

隨著計算機技術(shù)的快速發(fā)展和普及，有限元方法已迅速從工程結(jié)構(gòu)強度分析計算擴展到幾乎所有的科學技術(shù)領(lǐng)域，成為一種應(yīng)用廣泛并且實用高效的數(shù)值分析方法。早期有限元分析的研究重點在于推導新的高效率求解方法和高精度的單元。隨著數(shù)值分析方法的逐步完善和計算機運算速度的飛速提高，整個計算系統(tǒng)用于求解運算的時間越來越短，而數(shù)據(jù)準備和運算結(jié)果的表現(xiàn)問題卻日益突出[1]。網(wǎng)格剖分作為建立有限元模型的一個重要環(huán)節(jié)，要求考慮的問題多，需要的工作量大，不同的網(wǎng)格劃分方式會對計算規(guī)模、計算結(jié)果和計算精度產(chǎn)生很大的影響。故而，對有限元網(wǎng)格剖分的研究十分必要。

有限元分析的最終目的是要還原一個實際工程系統(tǒng)的數(shù)學行為特征，即分析必須是針對一個物理原型的準確的數(shù)學模型。從廣義上講，模型包括所有的節(jié)點、單元、材料屬性、幾何特性、初始條件、邊界條件等，以及其他用來表現(xiàn)這個物理系統(tǒng)的特征。從狹義上講，模型生成僅指用節(jié)點和單元表示空間體域以及實際系統(tǒng)連接的生成過程，即網(wǎng)格剖分。

在建立有限元模型過程中，不管從廣義還是從狹義上講，都涉及網(wǎng)格剖分問題。曾經(jīng)有人作過統(tǒng)計，在數(shù)值分析的三個階段中，前處理約占總時間的40%～60%，數(shù)值求解約占5%～20%，計算結(jié)果后處理約占30%[2]。如果純粹采用人工方法進行分析對象的離散化工作，勢必需要花費大量的時間，而且當模型復(fù)雜時，還容易出錯。前處理工作比較繁瑣，但是又十分重要，它是進行有限元正確分析的基礎(chǔ)。因此，開展更好的分析對象離散化網(wǎng)格劃分工作，對于數(shù)值分析工作者來講，具有非常重要的意義。

1 網(wǎng)格剖分要求

有限元網(wǎng)格生成就是將工作環(huán)境下的物體離散成單元的過程，常用的單元包括：一維桿單元及集中質(zhì)量元，二維三角形和四邊形單元，三維四面體、五面體、金字塔形、六面體單元。其邊界形狀主要有直線型、直面型、曲線型和曲面型。對于邊界為曲線(面)型的單元，有限元分析要求各邊或面上有若干節(jié)點，這樣既可保證單元的形狀，同時又能提高求解精度、準確性和加快收斂速度。不同維數(shù)的同一物體可劃分為多種單元混合而成的網(wǎng)格。

網(wǎng)格劃分應(yīng)該遵循以下原則。①合法性。一個單元的節(jié)點不能落入其他單元內(nèi)部，在單元邊界上的節(jié)點均應(yīng)作為單元的節(jié)點，不可丟棄。②相容性。單元必須落在待分區(qū)域內(nèi)部，不能落入外部，且單元并集等于待分區(qū)域。③協(xié)調(diào)性。單元上的力和力矩能夠通過節(jié)點傳遞給相鄰單元。為保證單元協(xié)調(diào)，必須滿足：一個單元的節(jié)點必須同時也是相鄰單元的節(jié)點，而不應(yīng)是內(nèi)點或邊界點；相鄰單元的共有節(jié)點具有相同的自由度性質(zhì)，即自由度必須匹配。④逼近精確性。待分區(qū)域的頂點(包括特殊點)必須是單元的節(jié)點，待分區(qū)域的邊界(包括特殊邊及面)被單元邊界所逼近。⑤良好的單元形狀[1]。單元最佳形狀是正多邊形或正多面體。⑥良好的劃分過渡性[1]。單元之間過渡應(yīng)相對平穩(wěn)，否則將影響計算結(jié)果的準確性甚至使有限元計算無法進行下去。⑦網(wǎng)格劃分的自適應(yīng)性[1]。在幾何尖角處、應(yīng)力、溫度等變化大的地方網(wǎng)格應(yīng)密，其他部位應(yīng)較稀疏，這樣可以保證計算結(jié)果精確可靠。⑧一致性。對于相連的兩個二次單元，單元角點只能與單元角點連接，而不能與相鄰單元的中間節(jié)點連接；相鄰單元的公共邊應(yīng)具有相同的節(jié)點數(shù)，當采用混合單元(線性單元與高階單元)類型時有必要從一個單元中除去中間節(jié)點。另外，在動力分析中，沖擊波傳播問題不推薦使用二次單元。

2 網(wǎng)格剖分準備與剖分方法

2.1 網(wǎng)格剖分準備工作

(1)確定合適的網(wǎng)格密度。在數(shù)值分析中經(jīng)常碰到的問題是：單元網(wǎng)格應(yīng)剖分得如何細致才能獲得合理的結(jié)果。對于此問題，可借助于以下一些技術(shù)解決：①利用自適應(yīng)網(wǎng)格剖分產(chǎn)生可以滿足某種誤差估計準則的網(wǎng)格。②與先前獨立得出的實驗分析結(jié)果或已知解析解進行對比。對已知結(jié)果和計算結(jié)果偏差過大的地方進行網(wǎng)格細化。③執(zhí)行一個認為是合理的網(wǎng)格剖分的初始分析過程，再在危險區(qū)域利用兩倍多的網(wǎng)格重新分析并比較兩者的結(jié)果。如果這兩者給出的結(jié)果幾乎相同，則網(wǎng)格是足夠的。如果產(chǎn)生了顯著不同的結(jié)果，應(yīng)該繼續(xù)細化網(wǎng)格直到隨后的剖分獲得了近似相等的結(jié)果。④如果細化網(wǎng)格測試顯示只有模型的一部分需要更細的網(wǎng)格，可以對模型使用子模型以放大危險區(qū)域。網(wǎng)格剖分密度很重要，如果網(wǎng)格過于粗糙，那么結(jié)果可能包含嚴重的錯誤，如果網(wǎng)格過于細致，將花費過多的計算時間，浪費計算資源，而且模型過大可能會導致不能在自己的計算機上運行。因而，在生成模型前應(yīng)仔細考慮網(wǎng)格密度問題。

(2)單元形狀與類型的選擇[3]。從單元幾何形狀看，一維分析有兩點或三點線單元，二維分析有三角形或四邊形單元，三維分析有四面體或六面體單元。四邊形單元可以退化成三角形單元，六面體單元可以退化成五面體(楔形或金字塔形)單元或四面體單元，這取決于所分析問題的維數(shù)。從數(shù)值積分方案看，有全積分和降階積分之分，取決于分析問題所要求的精度，全積分方案精度高，降階積分方案精度低。從分析問題的類型看，有結(jié)構(gòu)分析單元與非結(jié)構(gòu)分析單元，不同的節(jié)點自由度和不同的場問題控制方程需要采用不同的單元類型，具體有平面應(yīng)力、平面應(yīng)變、二維梁、軸對稱殼、軸對稱實體、三維殼、三維實體、三維梁、熱傳導單元等。從單元的階次看，有線性單元(無中間節(jié)點)和二次單元(帶中間節(jié)點)，利用二次單元分析的精度較高，但是由于節(jié)點數(shù)猛增，會增加計算成本，而且所需存儲空間也會成倍增加。

2.2 網(wǎng)格剖分方法[3]

(1)網(wǎng)格直接生成法(直接建模法)。網(wǎng)格直接生成技術(shù)并不僅僅只是手動添加單元，而是以單元作為基本構(gòu)造模塊。先用一個個比較大的單元組成一個很粗糙的(待細分)網(wǎng)格體，然后通過特定的工具對這個比較粗糙的網(wǎng)格再進行精細的重新劃分，達到所要求的精度。這個生成過程特別適合那些幾何模型簡單的問題，生成方法簡單，大體可分為三個步驟：①生成粗糙構(gòu)造模塊，并對內(nèi)部進行精細劃分，即生成坐標基本合適和具有完全合格的連通性的網(wǎng)格模型，然后在需要的地方重新劃分和細分那些單元；②使邊界坐標符合要求，即把邊界節(jié)點坐標更改正確；③重新調(diào)整內(nèi)部網(wǎng)格，即重新分配內(nèi)部的坐標來生成合理形狀或者“放松”的單元。

(2)由幾何實體生成網(wǎng)格法(實體建模法)。通過幾何實體生成網(wǎng)格，能夠?qū)⒂蓭缀卧孛枋龅奈锢砟Ｐ妥詣与x散成有限單元，應(yīng)用這種技術(shù)生成的模型的基本構(gòu)成模塊是幾何體而不是網(wǎng)格體。用來表征幾何體的幾何元素是點、線、面、體，幾何模型必須完全建立好之后才能被剖分成網(wǎng)格模型。由幾何模型生成網(wǎng)格的方法有很多好處，最主要的優(yōu)點在于復(fù)雜模型通過這種方法容易生成，而通過直接生成法往往很難處理，而且還容易出錯。這是該方法適用范圍更廣的重要原因。通過這種方法，有時還可以使用從其他CAD軟件傳輸過來的模型進行操作，可大大提高效率。

3 有限元網(wǎng)格剖分算法

網(wǎng)格剖分算法主要包含以下幾種：拓撲分解法、節(jié)點連元法、網(wǎng)絡(luò)模板法、映射法、幾何分解法、基于柵格法、空間編碼法[1,4-17]。目前，在許多商業(yè)軟件中，這些方法基本上是混合使用的，很少有單獨使用的，并且這些方法與現(xiàn)代計算機技術(shù)結(jié)合緊密。

3.1 拓撲分解法

拓撲分解法是由Wordenweber[15]提出的，用于求解二維平面問題，現(xiàn)已推廣至三維空間。該方法用一種三角化算法將目標用盡量少的三角形完全分割覆蓋，這些三角形主要是由目標的拓撲結(jié)構(gòu)決定，這樣目標的復(fù)雜拓撲結(jié)構(gòu)被分解成簡單的三角形拓撲結(jié)構(gòu)。該方法后來被發(fā)展為普遍使用的目標初始三角化算法，用來實現(xiàn)從實體表述到初始三角化表述的自動轉(zhuǎn)換。拓撲分解法原理簡單，引入的算子概念使程序易于實現(xiàn)模塊化，處理容易。但是該方法只從拓撲關(guān)系入手，不考慮幾何因素，因此難以保證網(wǎng)格質(zhì)量，而且檢測量很大，對包含曲面的三維形體也難以處理。

Woo等[16]提出了另外一種基于拓撲分解法的有限元網(wǎng)格自動生成算法，試圖解決三維實體的有限元網(wǎng)格生成問題，但其網(wǎng)格質(zhì)量難以得到保證。

3.2 節(jié)點連元法

節(jié)點連元法主要包括兩個步驟：節(jié)點生成和單元生成。首先在待分區(qū)域內(nèi)生成一定數(shù)目的節(jié)點，然后通過適當?shù)乃惴ㄟB接節(jié)點生成有限元單元。常用的算法有Delaunay三角剖分法(簡稱DT法)和推進波前法(advancing front technique，AFT)。

DT法是目前最流行的通用的全自動網(wǎng)格生成方法之一。其最大優(yōu)點是遵循“最小角最大”和“空球”準則?！白钚〗亲畲蟆笔侵冈诓怀霈F(xiàn)奇異性情況下，Delaunay三角剖分最小角之和均大于任何非Delaunay剖分所形成的三角形最小角之和。DT法自動避免了生成小內(nèi)角的長薄單元，特別適用于有限元網(wǎng)格生成。而“空球”準則是指在剖分的任意三角形單元或四面體單元的外接圓(二維)或外接球(三維)內(nèi)都不包含其他單元節(jié)點。DT法的計算效率與具體實現(xiàn)方法相關(guān)。

此網(wǎng)格劃分方法是先生成覆蓋區(qū)域的稀疏三角形單元，然后局部加密，再生成所需密度的三角形網(wǎng)格。所生成的單元形態(tài)趨向于等邊三角形。DT法充分考慮了幾何形狀中存在的微小幾何特征，并能在微小幾何特征處劃分較細的單元。在不需要密集網(wǎng)格處，采用稀疏單元，疏密網(wǎng)格的過渡十分平滑。

雖然DT法既適用于二維域也適用于三維域，但直接的DT法只適用于凸域，不適用于非凸域，因此后來又發(fā)展了多種非凸域的Delaunay剖分。

近年來，推進波前法也已經(jīng)成為目前最流行的通用的全自動網(wǎng)格生成方法之一。其基本原理是：設(shè)區(qū)域的有向離散外邊界集和邊界前沿點集已經(jīng)確定，按某種條件沿區(qū)域邊界向區(qū)域內(nèi)部扣除三角形(四面體)直到區(qū)域為空集。AFT的關(guān)鍵技術(shù)有兩個：區(qū)域的邊界離散與內(nèi)部節(jié)點合理生成，扣除三角形條件。而扣除三角形條件有多種：最短距離條件、最大角條件、最大形狀質(zhì)量條件、最小外接圓條件等。

AFT可以全自動地在平面或曲面上生成網(wǎng)格，用戶可控制生成單元的幾何分類：四邊形或三角形，或者四邊形和三角形混合網(wǎng)格。這種自動計算網(wǎng)格的方法是一次生成一個單元，從區(qū)域的邊界向內(nèi)部逐漸生成全域網(wǎng)格。由它生成的網(wǎng)格同樣有著很好的幾何尺寸和形狀且疏密過渡平滑。當網(wǎng)格疏密過渡較劇烈時，它也同樣能夠生成高質(zhì)量的網(wǎng)格。

3.3 網(wǎng)格模板法

網(wǎng)格模板法生成單元主要分兩步(以三維實體為例)：第一，將待剖分網(wǎng)格的實體用適當大小的立方體(樹根)完全包容，按照“一化八”的原則遞歸離散，通過對每個八分塊分類，形成IN和NIO八分塊的并集，稱為網(wǎng)格模板模型。第二，根據(jù)得到的網(wǎng)格模板模型，再進行網(wǎng)格劃分處理。

網(wǎng)格模板法的優(yōu)點是網(wǎng)格生成完全自動，網(wǎng)格剖分速度快，非常適用于自適應(yīng)網(wǎng)格生成；主要缺點是邊界單元形狀難于完全保證。另外，網(wǎng)格模板法對物體的方向較敏感。

3.4 映射法

映射法的基本思想是：通過適當?shù)挠成浜瘮?shù)將待剖分物理域映射到參數(shù)空間中形成規(guī)則參數(shù)域；對規(guī)則參數(shù)域進行網(wǎng)格剖分；將參數(shù)域的網(wǎng)格反向映射回物理空間，從而得到物理域的有限元網(wǎng)格。

當前許多商用有限元網(wǎng)格生成器都以該算法為理論基礎(chǔ)。其主要特征是采用了“調(diào)配函數(shù)”概念。產(chǎn)生的網(wǎng)格整齊劃一，非常規(guī)則，但同時這也是該法的缺點，尤其是當待劃分區(qū)域不規(guī)則時，如瓶頸形狀，得到的網(wǎng)格形狀是畸變的。另外，映射法是非全自動方法，必須通過人工交互方式，將剖分對象先剖分成具有簡單拓撲關(guān)系的子域。但映射法處理曲面問題很有效。

映射法的優(yōu)點是：算法簡單、速度快、單元質(zhì)量好、密度可控制，它既可生成結(jié)構(gòu)化網(wǎng)格又可生成非結(jié)構(gòu)化網(wǎng)格，既可生成四邊形單元網(wǎng)格又可生成六面體單元網(wǎng)格，可用于曲面網(wǎng)格生成，可與形狀優(yōu)化算法集成等。

映射法一般可直接處理單連通域問題，但對于復(fù)雜多連通域問題，需要首先用手工或自動方法將待剖分域分解成幾何形狀規(guī)則的可映射子區(qū)域，然后在每個子區(qū)域內(nèi)應(yīng)用映射法。然而在實踐中仍有幾個難點需要克服：① 如何自動地將復(fù)雜的不可映射的待剖分域分解成簡單的可映射的子區(qū)域；② 如何滿足某些物理問題中對網(wǎng)格疏密過渡的要求；③ 如何滿足子區(qū)域之間的網(wǎng)格相容性要求。

3.5 幾何分解法

在產(chǎn)生節(jié)點的同時，也確定了節(jié)點之間的連接關(guān)系的網(wǎng)格剖分方法稱為幾何分解法。這類算法基于問題域的拓撲幾何描述，通過從域中逐個移去單元而生成有限元網(wǎng)格。它較多地考慮了待分域的幾何特征，確保生成質(zhì)量較好的網(wǎng)格單元，通常有兩種方法：遞歸法和迭代法。幾何分解法可實現(xiàn)從實體幾何描述到初始網(wǎng)格生成之間的自動轉(zhuǎn)換，并允許網(wǎng)格密度變化，但只能通過邊界點的分布來控制網(wǎng)格規(guī)模，網(wǎng)格質(zhì)量不高，且很難實現(xiàn)局部自適應(yīng)加密。

幾何分解法的最大優(yōu)勢是離散時考慮了網(wǎng)格的形狀和大小，因此，所生成的網(wǎng)格單元形狀和分布比較好。但是，這種方法自動化程度比較低，也不利于復(fù)雜件的網(wǎng)格生成。

3.6 基于柵格法

基于柵格法(grid-based approach)的基本剖分流程如下：首先用一組不相交的尺寸相同或不同的柵格(cells)覆蓋在目標區(qū)域上面，保留完全或部分落在目標區(qū)域之內(nèi)的柵格，刪除完全落在目標區(qū)域之外的柵格；然后對與物體邊界相交的柵格進行調(diào)整、剪裁、再分解等操作，使其更準確地逼近目標區(qū)域；最后對內(nèi)部柵格和邊界柵格(特別是后者)進行柵格級的網(wǎng)格剖分，進而得到整個目標區(qū)域的有限元網(wǎng)格。

這種方法預(yù)先產(chǎn)生網(wǎng)格模板，然后將要進行網(wǎng)格化的物體加到其上，并在實體內(nèi)部盡可能多地填充規(guī)則的長方體或正方體網(wǎng)格，在實體的邊界上根據(jù)實體邊界的具體特征更改網(wǎng)格的形狀和相互連接關(guān)系，使邊界上的單元盡可能無限地逼近物體的邊界形狀。這種方法能實現(xiàn)網(wǎng)格生成的自動化，網(wǎng)格的生成速度也非?？?，能夠生成的單元類型很多，劃分簡單，效率較高。其最大缺點是物體邊界單元的質(zhì)量較差；另一個缺點是所生成的單元尺寸相近，網(wǎng)格密度很難得到控制。

柵格法首先用交互方式將物體劃分為形狀簡單的子區(qū)域，每個子區(qū)域分別用定形的網(wǎng)格模板作為規(guī)整的部分，再采取適當?shù)拇胧?，使得相鄰子區(qū)域在結(jié)合面共享公共的節(jié)點，并且網(wǎng)格相容。

3.7 空間編碼法

空間編碼法有兩個本質(zhì)屬性：階梯結(jié)構(gòu)和空間可訪問性。該法可以實現(xiàn)與實體造型系統(tǒng)的集成，并且容易精整網(wǎng)格質(zhì)量。大多數(shù)實體造型系統(tǒng)都采用樹形數(shù)據(jù)結(jié)構(gòu)進行幾何和拓撲描述，在此基礎(chǔ)上，Yerry等[18]提出了修正的八叉樹法等有限元網(wǎng)格生成算法。

基于修正的八叉樹法的空間編碼法在問題域內(nèi)部容易生成高質(zhì)量的單元，但是邊界單元需要進一步處理，以免所生成的單元因質(zhì)量太差而不適合有限元分析。也有學者將這種方法劃歸于基于柵格法。

4 網(wǎng)格形狀質(zhì)量判定及指標

4.1 網(wǎng)格質(zhì)量要求

如果單元都是理想的形狀(三角形單元等邊，四邊形都是正方形，六面體都是立方體等)，那么在計算單元剛度矩陣的時候誤差和錯誤就會很少出現(xiàn)，然而在整個模型中都用理想形狀的單元來離散幾乎是不可能的事情，尤其是模型非常復(fù)雜時，難以全用規(guī)則單元來剖分，因此必須在模型不同位置合理設(shè)置不同形狀的網(wǎng)格，改善網(wǎng)格質(zhì)量。

網(wǎng)格質(zhì)量對于數(shù)值分析的精度有十分重要的影響，特別對于具有復(fù)雜形狀的分析對象，尤為重要。對于復(fù)雜幾何實體和殼，首先要保證網(wǎng)格單元與幾何體嚴格對應(yīng)；其次保證單元具有高質(zhì)量。很多時候在進行數(shù)值計算過程中，會出現(xiàn)單元質(zhì)量不合格的警告信息，嚴重時會因出錯而導致計算中斷。因此，建立網(wǎng)格質(zhì)量標準體系，對于數(shù)值分析具有很重要的意義。

一般而言，對于不同的數(shù)值分析內(nèi)容網(wǎng)格質(zhì)量的標準有所不同。對于一般的熱傳導分析，可以適當降低網(wǎng)格質(zhì)量標準，而對于一些非線性問題，如大變形、接觸、瞬態(tài)高速沖擊等分析問題，需要高質(zhì)量的網(wǎng)格。

由于有限元網(wǎng)格的質(zhì)量直接影響到數(shù)值求解的精確度和正確性，自動生成的初始網(wǎng)格的質(zhì)量并不總是令人滿意的，所以通常要對初始網(wǎng)格的質(zhì)量進行改進或優(yōu)化。網(wǎng)格質(zhì)量的改進可分為兩個方面[18]：一是拓撲關(guān)系的調(diào)整；二是節(jié)點位置的調(diào)整。

檢查網(wǎng)格質(zhì)量，首先應(yīng)該檢查是否有重復(fù)的節(jié)點和單元。如果在同一位置需要建立兩個節(jié)點或單元，應(yīng)務(wù)必小心。在絕大部分商業(yè)工程分析軟件中，具有在一定公差范圍內(nèi)消除重復(fù)節(jié)點的功能，如果需要在同一位置產(chǎn)生不同節(jié)點(如接觸問題)，應(yīng)避免使用此功能。

4.2 一維網(wǎng)格質(zhì)量評價指標

一維網(wǎng)格質(zhì)量評價指標包括自由端點和剛性單元，即檢查網(wǎng)格內(nèi)是否存在自由端點和剛性單元。其中，自由端點主要檢查是否存在自由端點或自由節(jié)點(即與其他單元不相連)，在一維單元容易出現(xiàn)這個問題，如質(zhì)量集中單元等。剛性單元主要檢查是否具有形成環(huán)狀的剛性單元。

4.3 二維網(wǎng)格質(zhì)量評價指標

二維單元的幾何形狀主要是三角形和四邊形，主要質(zhì)量指標包括：單元長度、翹曲角、單元邊長比、內(nèi)角大小、扭曲角、雅可比比率(Jacobian ratio)等。

三角形單元主要檢查：單元長度、長寬比、扭曲角和內(nèi)角大小。

四邊形單元主要檢查：單元長度、翹曲度、長寬比、扭曲角、雅可比比率、弦偏離度。

各項技術(shù)質(zhì)量指標如下所述。

(1)單元邊長比γAR。γAR為單元最大邊長與最小邊長之比，適用于所有單元。對單元尺寸應(yīng)進行適當控制，既保證計算精度，又不浪費資源。理想的單元是單元邊長比為1，對線性單元來說，可接受的單元邊長比的范圍是0<γAR<3，對二次單元來說，可接受的單元邊長比范圍是0<γAR<10。對于同樣形狀的單元來說，高階單元對于邊長比沒有線性單元對于邊長比那么敏感。單元在非線性分析中對于邊長比的敏感程度要比在線性分析中對于邊長比的敏感程度高。如果一個問題在某一方向應(yīng)力梯度很大，單元有可能需要相當大的邊長比，最小邊放在梯度最大的地方。這是因為在一個單元內(nèi)，如果某一邊的梯度很大，這一邊又很長，那么誤差就會很大。

(2)三角形單元內(nèi)角。即三角形三個內(nèi)角大小。

(3)三角形單元扭曲角。這一指標表征了單元在單元面內(nèi)的扭曲程度。其定義為：對應(yīng)邊中點連線的夾角中最小角的余角，即三角形單元扭曲角θskew=90°-min(α1,α2,α3)，α1、α2、α3為中內(nèi)角，見圖1。另外還有一種定義：單元相鄰邊夾角與60°之間的差值。

圖1 三角形單元扭曲角定義

(4)四邊形單元扭曲角。該指標的定義為：對應(yīng)邊中點連線的夾角中最小角的余角，即四邊形單元扭曲角θskew=90°-min(δ1,δ2)，見圖2。另外一種定義是：單元相鄰邊夾角與90°之間的差值。

圖2 四邊形單元扭曲角定義

(5)四邊形單元翹曲角。該指標表征了單元在單元的面外的翹曲程度，面外翹曲發(fā)生在單元面的節(jié)點不共面的時候。其定義如下：依次沿對角線將四邊形分為兩個三角形，尋找這兩個三角形所在面構(gòu)成夾角的最大值，該角即為翹曲角，即θwarp=max(α1，α2)，見圖3。

圖3 四邊形單元翹曲角

(6)弦偏離度。即單元各邊中點與各點在對應(yīng)邊上的投影點的距離值，見圖4中的L1、L2。

圖4 四邊形單元弦偏離度

(7)雅可比比率。即單元內(nèi)各個積分點Jocabian行列式值中的最小值與最大值之比，見圖5。計算公式如下：

(1)

式中，JR為雅可比比率；|J|min、|J|max分別為最小和最大雅可比行列式值。

且

|J|(-1，-1)=(x2-x1)(y4-y1)-
(x4-x1)(y2-y1)=l1l4sinθ1

(2)

式中，x2、x1、y4、y1、x4、x1、y2、y1、l1、l4、θ1參見圖5c。

(a)規(guī)則四邊形單元積分點示意圖(b)任意直邊四邊形單元積分點示意圖

(c)任意直邊四邊形單元整體坐標

(d)任意直邊四邊形單元局部(自然)坐標圖5 四邊形單元雅可比計算示意圖

4.4 三維網(wǎng)格質(zhì)量評價指標

三維單元質(zhì)量檢查指標與二維單元質(zhì)量檢查指標大同小異，但有些指標不一樣，如在四面體中，邊長比定義為單元最長邊與最短高之間的比值，見圖6。對于六面體單元，單元質(zhì)量檢查指標與二維單元差不多，但對于四面體單元，還需要另外檢查如下幾個指標：

(1)四面體單元坍塌(collapse)值。如圖7所示，其計算公式如下：

Tcollapse=min(hi/sqrt(Ai))/1.24i=1,2,3,4

式中，hi為各個頂點到對應(yīng)面的距離值；Ai為對應(yīng)面的面積；sqrt(·)為取平方根運算的函數(shù)。

圖6 四面體單元邊長比示意圖圖7 四面體單元坍塌值計算示意圖

(2)四面體單元的體積扭曲(skew)值。對于任意一個四面體單元，定義一個過該四面體四個頂點的外接球體，如圖8所示，再依照球體的半徑，計算出一個理想四面體的體積，該體積假定為Videal，實際四面體單元的體積為Vactual，參照理想四面體的體積，按照下面公式計算，就可以得到四面體單元的扭曲值：

Tvolumetric_skew=(Videal-Vactual)/Videal

式中，r為外接球的半徑。

圖8 四面體單元體積扭曲值計算示意圖

5 有限元網(wǎng)格質(zhì)量優(yōu)化

5.1 網(wǎng)格優(yōu)化

為了保證有限元分析結(jié)果令用戶滿意，有限元網(wǎng)格通常應(yīng)具備以下條件[8]：①所有單元接近理想形狀；②主要變量(溫度、速度等)變化梯度較大的地方網(wǎng)格密度較大；③粗細網(wǎng)格之間過渡均勻。但通常情況下，在有限元網(wǎng)格自動生成器所產(chǎn)生的網(wǎng)格中總存在一些畸形單元，問題域越復(fù)雜，畸形網(wǎng)格所占比例越大。

網(wǎng)格優(yōu)化目的是改變網(wǎng)格質(zhì)量，提高計算精度。目前主要存在兩種優(yōu)化方法[8,19]：①網(wǎng)格精整。根據(jù)網(wǎng)格密度和計算結(jié)果的要求對網(wǎng)格進行細分；②光滑處理。保持網(wǎng)格拓撲關(guān)系不變，通過攝動單元節(jié)點位置來改善網(wǎng)格質(zhì)量。

Laplacian光滑處理技巧是應(yīng)用得最早，也是最成熟的一種優(yōu)化方法[20]。其核心內(nèi)容為：保持網(wǎng)格拓撲關(guān)系不變，將整個內(nèi)部節(jié)點的位置攝動到由其相鄰節(jié)點組成的多邊形的質(zhì)心處，使每個單元更接近于理想形狀。將這個攝動過程遍歷所有內(nèi)部節(jié)點若干次，可較大幅度地提高網(wǎng)格質(zhì)量。經(jīng)過Laplacian光滑處理的網(wǎng)格，不再具有Delaunay三角剖分的基本性質(zhì)。

近年來出現(xiàn)了各種各樣的網(wǎng)格優(yōu)化方法，自適應(yīng)網(wǎng)格精整方法即是其中之一。首先用較稀疏的網(wǎng)格進行有限元分析，根據(jù)計算結(jié)果，通過某種誤差指示器判斷哪些區(qū)域的網(wǎng)格需要精整，經(jīng)局部自適應(yīng)網(wǎng)格精整后再對該區(qū)域進行有限元分析。這樣的局部網(wǎng)格精整過程可反復(fù)進行多次，直到計算結(jié)果滿足用戶要求為止。

5.2 自適應(yīng)網(wǎng)格方法

在目前很多數(shù)值計算中都需要用到網(wǎng)格自適應(yīng)技術(shù)，如激波產(chǎn)生的高速飛行、材料加工等大變形問題、瞬態(tài)高速沖擊、炸藥爆轟聚能射流及侵徹等。在有限元分析中，網(wǎng)格自適應(yīng)是在現(xiàn)有的網(wǎng)格基礎(chǔ)上，根據(jù)有限元計算結(jié)果估計計算誤差、重新剖分網(wǎng)格和再計算的循環(huán)過程。當計算誤差達到預(yù)定值時，自適應(yīng)過程結(jié)束。因此，有效的誤差估計和良好的自適應(yīng)網(wǎng)格生成是自適應(yīng)有限元分析兩大關(guān)鍵技術(shù)[9,21-23]。就目前國內(nèi)外研究來看，自適應(yīng)網(wǎng)格生成從大的方面可分為網(wǎng)格增加技術(shù)和網(wǎng)格重劃分技術(shù)兩類。

網(wǎng)格增加技術(shù)主要依靠增加自由度總數(shù)來提高有限元分析精度。目前，主要采用三種類型方法：h型、 p型和 h-p型[22]。h型使用特別廣泛，網(wǎng)格模板模型的網(wǎng)格改進正是利用該法。p型是在保持網(wǎng)格劃分不變的情況下，通過提高插值函數(shù)的階數(shù)獲得更高的求解精度。h-p型是將h型和p型結(jié)合的一種方法。h-p型雖然實現(xiàn)不容易，但它卻可使收斂速率明顯加快。實踐表明，在獲得同一精度時，上述三種類型收斂速率是按h型、p型、h-p型順序增大的。

網(wǎng)格重劃分技術(shù)是根據(jù)現(xiàn)有的網(wǎng)格并配合誤差估計確定新的節(jié)點分布，然后重新劃分網(wǎng)格，再計算，重復(fù)上述過程直到求解精度達到預(yù)定目標為止的過程[21]。目前，網(wǎng)格重劃分技術(shù)在平面區(qū)域已得到了較好的實現(xiàn)。從理論上講，該原理可擴充到三維實體域，但由于三維實體域難以完全自動地用等節(jié)點密度曲面來分割任意實體，因此在三維域的擴充至今仍未實現(xiàn)。實踐表明，網(wǎng)格重劃分技術(shù)比網(wǎng)格增加技術(shù)具有更多的優(yōu)點，如收斂速度快、網(wǎng)格單元形狀穩(wěn)定等。

6 網(wǎng)格剖分前處理軟件及其應(yīng)用狀況[24]

有限元技術(shù)的特點決定其前處理過程為一個相對獨立而又十分重要的部分。隨著眾多有限元商業(yè)軟件的出現(xiàn)，前處理技術(shù)也得到了迅猛的發(fā)展，同時很多專業(yè)的前處理軟件也應(yīng)運而生。目前，在國際上被認可的前處理軟件主要包括Altair公司的HyperMesh，MSC公司的Patran，EDS公司的FEMAP，SamTech公司的Samcef/Field，CAE-Beta公司的ANSA，CFDRC公司的CFD-GEOM和CFD-MicroMesh，TrueGrid，F(xiàn)luent軟件中的Gambit等，還有一些大型有限元商業(yè)軟件自帶有網(wǎng)格剖分器，如Marc/Mentat、Ansys/PreProcesser等。在一般情況下，前處理軟件都與CAD軟件具有良好的接口，且與眾多的有限元求解器結(jié)合，可使用戶更快、更方便地解決問題。一些大型企業(yè)都采用了適應(yīng)自己需求的前處理軟件。

當前，應(yīng)用最廣泛的前處理軟件首推HyperMesh，它是一款高效率的有限元前處理軟件，可與大多數(shù)的有限元分析軟件搭配使用，如Marc、Nastran、Abaqus、Ansys、Ls-Dyna等。HyperMesh主要用于汽車行業(yè)，它已經(jīng)成為全球汽車行業(yè)的標準配置，幾乎所有的整車廠商都在使用。同時HyperMesh也廣泛進入各行各業(yè)，如航空、航天、通用機械與日用品等行業(yè)。

MSC公司的Patran軟件是一個集成的并行框架式有限元前處理及分析仿真系統(tǒng)，最早由美國宇航局(NASA)倡導開發(fā)，是工業(yè)領(lǐng)域最著名的軟件系統(tǒng)。其開放式、多功能的體系結(jié)構(gòu)可集工程分析、結(jié)果評估、用戶化設(shè)計和交互圖形界面于一身，構(gòu)成一個完整的CAE集成環(huán)境。Patran的用戶主要集中在航空、航天、汽車和通用機械等領(lǐng)域。

FEMAP是一個純Windows風格的、非常易于使用的高性能有限元前處理軟件。FEMAP提供給工程師和分析人員一個可以容易、精確、有效地操控復(fù)雜模型的前處理手段。FEMAP從高級梁造型、中面提取和高級網(wǎng)格劃分，到功能卓越的直接訪問CAD工具和簡化工具，都提供了有效的方法。

Samcef/Field系列軟件的前處理工具是一個獨立的圖形環(huán)境，具有幾何建模、讀取主流CAD模型和驅(qū)動Samcef線性求解器的能力。

ANSA(automatic net-generation for structure analysis)是希臘BETA CAE System公司的軟件產(chǎn)品，是世界上廣泛應(yīng)用的、功能強大的前處理系統(tǒng)和三維網(wǎng)格軟件工具。ANSA起源于汽車工業(yè)領(lǐng)域，主要是為了滿足有限元前處理對時間的需求。ANSA已成為工業(yè)標準，而且在汽車、航空航天和化工等工業(yè)領(lǐng)域應(yīng)用廣泛。

CFD-MicroMesh是為微電子和微機電工業(yè)領(lǐng)域的特殊需求而開發(fā)的軟件系統(tǒng)。它直接從EDA布局圖和初始設(shè)計開始，是自動化的三維幾何創(chuàng)建和網(wǎng)格生成工具。

TrueGrid是美國XYZ Scientific Applications公司推出的專業(yè)通用的網(wǎng)格劃分前處理軟件，支持大部分有限元分析及計算流體動力學(CFD)軟件，可以方便快速生成優(yōu)化的、高質(zhì)量的、多塊結(jié)構(gòu)的六面體網(wǎng)格模型。TrueGrid支持一般CAD所輸出的幾何形狀，如AutoCAD、Pro/E、I-Deas等；支持多款當今主流的分析軟件，如Ansys、Abaqus、Adina、Ls-Dyna、Autodyn、Marc、Nastran和流體動力學分析軟件，如與Fluent、AutoReagas、CFX、CFdesign、Star-CD、Phoenics、Numeca、Tascflow等具有接口。

Gambit是幫助分析者和設(shè)計者建立并網(wǎng)格化計算流體力學模型和其他科學應(yīng)用而設(shè)計的一個軟件包，是面向CFD分析的高質(zhì)量的前處理器，其主要功能包括幾何建模和網(wǎng)格生成。由于Gambit所具有的強大功能，在目前所有的CFD前處理軟件中，Gambit穩(wěn)居上游。Gambit通過它的用戶界面(GUI)來接受用戶的輸入，能直接建立模型、網(wǎng)格化模型、指定模型區(qū)域大小等，這對很多的模型應(yīng)用已足夠。

目前，國外很多廠商主要以HyperMesh和ANSA作為劃分網(wǎng)格與前處理的工具標準。

7 網(wǎng)格剖分實際案例

利用HyperMesh軟件平臺，筆者進行了二次開發(fā)，對幾個復(fù)雜模型進行了網(wǎng)格剖分。該軟件具有先進的網(wǎng)格剖分與網(wǎng)格優(yōu)化算法，能自動對網(wǎng)格進行光滑、優(yōu)化處理，改善網(wǎng)格質(zhì)量，能進行完備的網(wǎng)格質(zhì)量檢查，具備強大的網(wǎng)格編輯功能，可以方便調(diào)整、編輯網(wǎng)格，以改善網(wǎng)格質(zhì)量。

實例一汽車發(fā)動機活塞的網(wǎng)格剖分(圖9),使用了網(wǎng)格生成映射法和直接建模方法。公司提供如下的網(wǎng)格質(zhì)量與技術(shù)要求：①總的單元數(shù)量為50 000～60 000個；②在缸體內(nèi)壁面上不能出現(xiàn)三角形單元；③90%以上的單元是八節(jié)點六面體單元；④單元必須是六面體單元和五面體楔形單元；⑤網(wǎng)格質(zhì)量指標為θwarp<30°，γAR<5，θskew<60°，JR>0.3；⑥不能出現(xiàn)破裂和T形連接。

圖9 汽車發(fā)動機汽缸活塞六面體網(wǎng)格模型

圖10 復(fù)雜模型六面體網(wǎng)格模型

實例二某復(fù)雜模型網(wǎng)格剖分(圖10)，綜合使用了幾何分解法、推進波前法和直接建模法。網(wǎng)格剖分要求為：①總的單元數(shù)量為20 000～30 000個；②95% 以上的單元是八節(jié)點六面體單元；③單元形狀必須是六面體單元或五面體楔形單元；④網(wǎng)格質(zhì)量指標為θwarp<40°，γAR<4，θskew<50°，JR>0.5；⑤不能出現(xiàn)破裂和T形連接。

實例三汽車轉(zhuǎn)向盤網(wǎng)格剖分(圖11)，使用DT法。其網(wǎng)格質(zhì)量與技術(shù)要求為：①全部用四面體單元；②單元數(shù)量為30 000～40 000個；③網(wǎng)格質(zhì)量指標為單元邊長小于6mm，γAR<4，θskew<50°，JR>0.5；④不能出現(xiàn)T形連接。

圖11 汽車方向盤四面體網(wǎng)格模型

8 網(wǎng)格剖分前處理技術(shù)發(fā)展趨勢

(1)通用算法的數(shù)據(jù)結(jié)構(gòu)與多種算法的聯(lián)合應(yīng)用[1]。在通用算法的研究方面，應(yīng)注意數(shù)據(jù)結(jié)構(gòu)的研究和多種算法的聯(lián)合應(yīng)用，提高核心算法的可靠性和幾何適應(yīng)性，達到速度與質(zhì)量之間的平衡，實現(xiàn)核心算法的黑箱化。

(2)根據(jù)所分析問題的物理特性對網(wǎng)格進行修改，從這個角度出發(fā)有兩個主要的問題：生成具有所有性質(zhì)或部分性質(zhì)的網(wǎng)格，對現(xiàn)有的網(wǎng)格進行修改以獲得所期望的性質(zhì)，其中自適應(yīng)有限元分析是一個最為主要的研究領(lǐng)域[25]。

(3)自適應(yīng)網(wǎng)格剖分技術(shù)。自適應(yīng)算法中通過自動加密網(wǎng)格來提高計算精度，這種算法有效而且收斂速度快，在各個領(lǐng)域內(nèi)自適應(yīng)算法將會有較大的發(fā)展，其主要問題在于如何進行誤差估算。這對于各個領(lǐng)域中的不同問題是不一樣的，這也是目前有限元前處理技術(shù)發(fā)展的主要方向之一。

(4)六面體網(wǎng)格自動剖分技術(shù)。幾十年來，眾多學者致力于六面體單元網(wǎng)格自動生成方法研究，但復(fù)雜三維實體的全六面體單元網(wǎng)格全自動生成問題始終未能獲得真正意義上的解決。近幾年來，全六面體網(wǎng)格自動生成再度成為焦點問題[26-33]。

(5)網(wǎng)格生成算法的并行化和分布化。并行化計算環(huán)境對于大規(guī)模、超大規(guī)模科學計算以及高端工程應(yīng)用是必需的，而分布式計算環(huán)境可作為一種中端工程應(yīng)用解決方案?，F(xiàn)有網(wǎng)格生成并行化或分布化算法在計算效率、內(nèi)存管理、生成單元質(zhì)量等方面還不夠完善，還有許多潛力可挖。另外，并行計算環(huán)境與分布式計算環(huán)境的控制軟件日趨成熟，這為算法的并行化、分布化開發(fā)提供了更強有力的技術(shù)保障。

9 結(jié)語

數(shù)值分析已經(jīng)與理論研究、實(試)驗研究成為三大研究技術(shù)手段之一，在計算機技術(shù)高度發(fā)展的今天，這種手段將會發(fā)揮越來越大的作用。采用離散化方法對研究對象進行網(wǎng)格剖分是數(shù)值模擬前處理的核心問題。研究網(wǎng)格剖分方法、建立網(wǎng)格質(zhì)量標準及技術(shù)要求體系對于有限元分析具有重要意義。在數(shù)值模擬應(yīng)用日益廣泛的今天，必須更加深入研究網(wǎng)格剖分技術(shù)方法，建立網(wǎng)格質(zhì)量標準及技術(shù)要求體系。

[1] 王明強，朱永梅，劉文欣.有限元網(wǎng)格劃分方法應(yīng)用研究[J].機械設(shè)計與制造，2004(1): 22-24.

[2] Shephard M S. Approaches to the Automatic Generation and Control of Finite Element Meshes[J]. Applied Mechanics Review,1998,4(4):169-185.

[3] 鄭巖，顧松東，吳斌.Marc 2001從入門到精通[M].北京：中國水利水電出版社，2003.

[4] 呂軍，王忠金，王仲仁.有限元六面體網(wǎng)格的典型生成方法及發(fā)展趨勢[J].哈爾濱工業(yè)大學學報，2001，33(4)：485-490.

[5] 古成中，吳新躍，有限元網(wǎng)格劃分及發(fā)展趨勢[J].計算機科學與探索，2008，2(3), 248-259.

[6] 關(guān)振群，宋超，顧元憲，等.有限元網(wǎng)格生成方法研究的新進展[J].計算機輔助設(shè)計與圖形學學報，2003，15(1)：1-14.

[7] Ho-Le K. Finite Element Mesh Generation Aethods:a Review and Classification[J].Computer-Aided Design,1988,20(1):27-38.

[8] 李俊，游理華.有限元網(wǎng)格自動生成算法的研究進展[J].機械研究與應(yīng)用，1998，11(4)：25-27.

[9] 李衛(wèi)民.有限元網(wǎng)格生成算法的評述[J].泰州職業(yè)技術(shù)學院學報，2004，4(1)：12-14.

[10] 胡恩球，張新訪，向文，等.有限元網(wǎng)格生成方法發(fā)展綜述[J].計算機輔助設(shè)計與圖形學學報，1997，9(4)：378-383.

[11] 張建華，葉尚輝.有限元網(wǎng)格自動生成典型方法與發(fā)展方向[J].技術(shù)研究，1996(2):28-31.

[12] 閔衛(wèi)東，唐澤圣.有限元網(wǎng)格劃分技術(shù)[J].計算機研究與發(fā)展，1995，32(7)：37-42.

[13] George P L.Automatic Mesh Generation. Applications to Finite Element Methods[M]. New York: Willey,1991.

[14] Lohner R, Juan R C. Parallel Advancing Front Grid Generation[C]//Proceedings of the 8th International Meshing Roundtable. Sout Lake Tahoe,CA,1999:67-74.

[15] Wordenweber B. Finite Element Mesh Generation[J]. Computer-Aided Design,1984,16(5)：15-20.

[16] Woo T C, Thomasma T. An Algorithm for Generating Solid Element in Objects with Holes[J]. Comp. Struct., 1984, 18(2)：30-36.

[17] Ruppert J, Seidel R. On the Difficulty of Triangulating Three-dimensional Nonconvex Polyhedra[J]. Discrete & Computational Geometry, 1992,7(3):227-253.

[18] Yerry M A,Shephard M S.Automatic Three Dimensional Mesh Generation by the Modified-octree Technique[J].Int. J. Num. Methods. Eng.，1984,20(11)：1965-1990.

[19] 黃志超，包忠詡，周天瑞，等.有限元網(wǎng)格劃分技術(shù)研究[J].南昌大學學報，2001，23(4)：25-32.

[20] Huang C Y. Recent Progress in Multiblock Hybrid Structured and Unstructured Mesh Generation[J].Computer Methods in Applied Mechanics and Engineering, 1997,150(1/4):1-24.

[21] Almeida R C, Feijo R A, Galeao A C,et al.Adaptive Finite Element Computational Fluid Dynamics Using an Anisotropic Error Estimator[J]. Computer Methods in Applied Mechanics and Engineering, 2000,182(4):379-400.

[22] Fuenmayor F J,Restrepo J L,Tarancn J E, et al. Error Estimation and H-adaptivere Finement in the Analysis of Natural Frequencies[J]. Finite Elements in Analysis and Design, 2001, 38(2):137-153.

[23] 單菊林.自適應(yīng)有限元網(wǎng)格生成算法研究與應(yīng)用[D].大連：大連理工大學，2007.

[24] 于開平，周傳月，譚惠豐,等.HyperMesh從入門到精通[M].北京：科學出版社，2005.

[25] 羅特軍，羅季軍，汪榴.有限元網(wǎng)格優(yōu)化方法[J].四川聯(lián)合大學學報，1999，3(3)：65-72.

[26] Chrisochoides N, Nave D. Simultaneous Mesh Generation and Partitioning for Delaunay Meshes[J]. Mathematics and Computers in Simulation,2000,54(4/5):321-339.

[27] Liu Shangsheng, Rajit G. Automatic Hexahedral Mesh Generation by Recursive Convexand Swept Volume Decomposition[C]//Proceedings of 6th International Meshing Roundtable.Washington: Sandia National Laboratories,1997:217-2311.

[28] Blacker T.Automated Conformal Hexahedral Meshing Constraints,Challenges and Opportunities[J]. Engineering with Computers, 2001, 17(3):201-210.

[29] Schneiders R, Bunten R. Automatic Generation of Hexahedral Finite Element Meshes[J].Computer Aided Geometric Design,1995,12:693-70.

[30] Schneiders R,Weiler F.Octree-based Generation of Hexahedral Element Mesher[C]//5th International Meshing Roundtable.Washington:Sandia National Laboratories,1996:205-216.

[31] Mller-Hannemann M. Quadrilateral Surface Meshes Without Self-intersecting Dual Cycles for Hexahedral Mesh Generation[J]. Computional Geometry, 2002,22(1/3):75-97.

[32] Takeo Taniguchi, Tomoaki Goda, Harald Kasper, et al. Hexahedra Mesh Generation of Complex Composite Domain[C]//Proceedings of the 5th International Conference on Grid Generation in Computationa Field Simulations. Mississippi:Mississippi State University, 1996:699-707.

[33] 關(guān)振群，單菊林，顧元憲.基于轉(zhuǎn)換模板的三維實體全六面體網(wǎng)格生成方法[J].計算力學學報，2005，22(1)：32-37.