有限元網格剖分與網格質量判定指標

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

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

0 引言

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

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

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

1 網格剖分要求

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

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

2 網格剖分準備與剖分方法

2.1 網格剖分準備工作

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

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

2.2 網格剖分方法[3]

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

(2)由幾何實體生成網格法(實體建模法)。通過幾何實體生成網格，能夠將由幾何元素描述的物理模型自動離散成有限單元，應用這種技術生成的模型的基本構成模塊是幾何體而不是網格體。用來表征幾何體的幾何元素是點、線、面、體，幾何模型必須完全建立好之后才能被剖分成網格模型。由幾何模型生成網格的方法有很多好處，最主要的優(yōu)點在于復雜模型通過這種方法容易生成，而通過直接生成法往往很難處理，而且還容易出錯。這是該方法適用范圍更廣的重要原因。通過這種方法，有時還可以使用從其他CAD軟件傳輸過來的模型進行操作，可大大提高效率。

3 有限元網格剖分算法

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

3.1 拓撲分解法

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

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

3.2 節(jié)點連元法

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

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

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

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

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

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

3.3 網格模板法

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

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

3.4 映射法

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

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

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

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

3.5 幾何分解法

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

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

3.6 基于柵格法

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

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

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

3.7 空間編碼法

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

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

4 網格形狀質量判定及指標

4.1 網格質量要求

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

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

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

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

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

4.2 一維網格質量評價指標

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

4.3 二維網格質量評價指標

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

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

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

各項技術質量指標如下所述。

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

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

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

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

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

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

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

圖3 四邊形單元翹曲角

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

圖4 四邊形單元弦偏離度

(7)雅可比比率。即單元內各個積分點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 三維網格質量評價指標

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

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

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

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

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

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

Tvolumetric_skew=(Videal-Vactual)/Videal

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

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

5 有限元網格質量優(yōu)化

5.1 網格優(yōu)化

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

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

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

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

5.2 自適應網格方法

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

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

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

6 網格剖分前處理軟件及其應用狀況[24]

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

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

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

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

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

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

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

TrueGrid是美國XYZ Scientific Applications公司推出的專業(yè)通用的網格劃分前處理軟件，支持大部分有限元分析及計算流體動力學(CFD)軟件，可以方便快速生成優(yōu)化的、高質量的、多塊結構的六面體網格模型。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是幫助分析者和設計者建立并網格化計算流體力學模型和其他科學應用而設計的一個軟件包，是面向CFD分析的高質量的前處理器，其主要功能包括幾何建模和網格生成。由于Gambit所具有的強大功能，在目前所有的CFD前處理軟件中，Gambit穩(wěn)居上游。Gambit通過它的用戶界面(GUI)來接受用戶的輸入，能直接建立模型、網格化模型、指定模型區(qū)域大小等，這對很多的模型應用已足夠。

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

7 網格剖分實際案例

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

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

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

圖10 復雜模型六面體網格模型

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

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

圖11 汽車方向盤四面體網格模型

8 網格剖分前處理技術發(fā)展趨勢

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

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

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

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

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

9 結語

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

[1] 王明強，朱永梅，劉文欣.有限元網格劃分方法應用研究[J].機械設計與制造，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] 呂軍，王忠金，王仲仁.有限元六面體網格的典型生成方法及發(fā)展趨勢[J].哈爾濱工業(yè)大學學報，2001，33(4)：485-490.

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

[6] 關振群，宋超，顧元憲，等.有限元網格生成方法研究的新進展[J].計算機輔助設計與圖形學學報，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] 李俊，游理華.有限元網格自動生成算法的研究進展[J].機械研究與應用，1998，11(4)：25-27.

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

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

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

[12] 閔衛(wèi)東，唐澤圣.有限元網格劃分技術[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] 黃志超，包忠詡，周天瑞，等.有限元網格劃分技術研究[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] 單菊林.自適應有限元網格生成算法研究與應用[D].大連：大連理工大學，2007.

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

[25] 羅特軍，羅季軍，汪榴.有限元網格優(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] 關振群，單菊林，顧元憲.基于轉換模板的三維實體全六面體網格生成方法[J].計算力學學報，2005，22(1)：32-37.