六面體單元生成方法及相關(guān)技術(shù)

1、3.1 引言有限元法是求解偏微分方程描述的連續(xù)體問(wèn)題的一種近似工程方法。為克服實(shí)際連續(xù)體問(wèn)題難于處理，它將分析區(qū)域離散化，將偏微分方程轉(zhuǎn)化為線性方，采用數(shù)值計(jì)算方法求出連續(xù)體問(wèn)題的近似解。用有限元法進(jìn)行工程分析的主要過(guò)程包括三個(gè)階段：（1）有限元模型的建立和數(shù)據(jù)準(zhǔn)備；（2）用分析計(jì)算；（3）分析結(jié)果的判斷和評(píng)定。迅速而合理地劃分有限元網(wǎng)格是完成有限元分析的前提和保證。有限元網(wǎng)格生成技術(shù)發(fā)展到現(xiàn)在，已經(jīng)出現(xiàn)了大量的、不同的實(shí)現(xiàn)方法 22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39。網(wǎng)格生成方法一般可以分為兩類：結(jié)構(gòu)化網(wǎng)格生成方法和非結(jié)構(gòu)網(wǎng)

2、格生成方法。所謂結(jié)構(gòu)化網(wǎng)格，嚴(yán)格的講，是指網(wǎng)格節(jié)點(diǎn)具有相同的鄰接單元的網(wǎng)格，即網(wǎng)格節(jié)點(diǎn)具有相同的拓?fù)浣Y(jié)構(gòu)。結(jié)構(gòu)化網(wǎng)格生成通常使用復(fù)合的迭代平滑技術(shù)，用邊界或物理區(qū)域來(lái)排列單元，劃分區(qū)域的邊界不能太復(fù)雜，以便能將區(qū)域分解為拓?fù)浣Y(jié)構(gòu)的塊。所謂非結(jié)構(gòu)化網(wǎng)格，放松了對(duì)節(jié)點(diǎn)的連接要求，允許一個(gè)節(jié)點(diǎn)屬于任意數(shù)量的單元。雖然四邊形和六面體單元也可以是非結(jié)構(gòu)網(wǎng)格，但通常非結(jié)構(gòu)網(wǎng)格是指三角形和四面體網(wǎng)格。目前二維區(qū)域網(wǎng)格自動(dòng)劃分發(fā)展較為成熟，許多商業(yè)都提供了成工具。而三維區(qū)域的網(wǎng)格自動(dòng)劃分，除個(gè)別商業(yè)提供了四面體網(wǎng)格的自動(dòng)劃分外，還不成熟，特別是形狀較規(guī)范的六面體網(wǎng)格的自動(dòng)生成。一般來(lái)說(shuō)，由程序自動(dòng)劃分的網(wǎng)格總

3、有一部分的單元的形狀不是很好，需要對(duì)網(wǎng)格的質(zhì)量進(jìn)行優(yōu)化。節(jié)點(diǎn)序號(hào)的的標(biāo)注方法直接影響剛度矩陣容量的大小從而影響計(jì)算的效率，因而需要對(duì)網(wǎng)格的節(jié)點(diǎn)重新，以得到具有較小帶寬的網(wǎng)格系統(tǒng)。本章主要了三維六面體網(wǎng)格生成的方法及其相關(guān)技術(shù)。3.2 結(jié)構(gòu)化網(wǎng)格生成結(jié)構(gòu)化網(wǎng)格生成方法主要分為兩類：代數(shù)法和 PDE（Partial Differential Equation，偏微分方法和保角法；PDE 方法包括橢圓系統(tǒng)方程、Poisson程）方法。代數(shù)法包括超限插值法、等參系統(tǒng)方程、拋物線系統(tǒng)方程和雙曲線系統(tǒng)方程方法。應(yīng)用最多的是等參法。結(jié)構(gòu)化網(wǎng)格的優(yōu)點(diǎn)是易于實(shí)現(xiàn)，在每個(gè)子區(qū)域內(nèi)網(wǎng)格可以得到很好的控制，生成規(guī)則的

4、結(jié)構(gòu)化網(wǎng)格，而且能生成曲面網(wǎng)格。缺點(diǎn)是對(duì)不規(guī)則的形體有時(shí)生成的網(wǎng)格質(zhì)量很差，需要事先根據(jù)所要產(chǎn)生的網(wǎng)格類型將目標(biāo)區(qū)域分割為一系列可的子區(qū)域，這一步通常需要手工完成，因而自動(dòng)化程度低，不適合網(wǎng)格全自動(dòng)生成。提供結(jié)構(gòu)化網(wǎng)格生成的有553DMAGGS（NASA，Langley & Lockhead，四邊形和六面體單元）、CAGI（ERC，Miss. Se，四邊形和六面體單元）、CSCMDO（CSC/NASA LaRC，四邊形單元）、EMC2（INRIA，四邊形單元）、GENIE+（ERC，Miss. Se，四邊形和六19面體單元）、ICEM CFD（ICEM CFD Engineering，四邊形和

5、六面體單元）、QulkGrid（Edge，四邊形單元）、TrueGrid（XYZ Scientific Applications，Inc.，四邊形和六面體單元）、VGM（NASA，Langley & Lockhead，四邊形和六面體單元）等。3.2.1 代數(shù)法網(wǎng)格生成代數(shù)法網(wǎng)格生成技術(shù)40是通過(guò)插值函數(shù)將理想的直角坐標(biāo)系表示的計(jì)算區(qū)域變換為實(shí)際物理區(qū)域來(lái)實(shí)現(xiàn)的（圖3-1）。它的具體做法是將計(jì)算區(qū)域的直角坐標(biāo)用均勻的間隔劃分為計(jì)算區(qū)域的網(wǎng)格，通過(guò)變換，計(jì)算區(qū)域上均勻分布的坐標(biāo)被為物理區(qū)域上的坐標(biāo)（圖3-2）。網(wǎng)格點(diǎn)數(shù)目的多寡并不影響變換的特性。從計(jì)算區(qū)域到物理區(qū)域的變換是通過(guò)一個(gè)向量函數(shù)完成的物

6、理區(qū)域計(jì)算區(qū)域zyx圖3-1計(jì)算區(qū)域和物理區(qū)域Fig. 3-1 Compuionaland physical計(jì)算網(wǎng)格物理網(wǎng)格x = x(, , )y = y(, , )z = z(, , )zyx圖3-2計(jì)算區(qū)域網(wǎng)格和物理區(qū)域網(wǎng)格Fig. 3-2 Gridshe compu ionalandhe physical20 x , , F , , y , , 0 1, 0 1, 0 1(3-1) z , , 多變量的坐標(biāo)變換式(3-1)使用無(wú)限插值(Transfiniteolation)來(lái)完成。在三維情況下，使用混合函數(shù)和與之相聯(lián)系的參數(shù)（特定點(diǎn)的位置及偏導(dǎo)數(shù)）來(lái)顯式地決定式(3-1)，然后通過(guò)對(duì)每

7、個(gè)單變量的循環(huán)完成無(wú)限插值。一般來(lái)說(shuō)，這些混合函數(shù)和參數(shù)都選取為區(qū)域邊界處的函數(shù)和參數(shù)。文獻(xiàn)40給出了幾種無(wú)限插值函數(shù)及其參數(shù)的選取。代數(shù)法無(wú)限插值網(wǎng)格生成法的特點(diǎn)是計(jì)算簡(jiǎn)單，可以采用中間變量的方法方便地控制網(wǎng)格的密度，對(duì)邊界簡(jiǎn)單的區(qū)域，可以生成質(zhì)量較高的網(wǎng)格，但缺點(diǎn)是不適應(yīng)復(fù)雜邊界的劃分，邊界不規(guī)則時(shí)生成的網(wǎng)格的質(zhì)量很差，并可能產(chǎn)生奇異性?？梢酝ㄟ^(guò)將劃分區(qū)域分解為子區(qū)域的方法，在子區(qū)域上應(yīng)用，可以在一定程度上克服這些缺點(diǎn)，但不易實(shí)現(xiàn)自動(dòng)劃分。3.2.2 偏微分方程法類似于代數(shù)法，偏微分方程法也是通過(guò)將曲的邊界為網(wǎng)格坐標(biāo)線，只不過(guò)這里使用的函數(shù)滿足某一類微分方程。以二維橢圓偏微分方程為例41，

8、滿足最大條件的橢圓型微分方可以為, Q , (3-2)xxyy右端的強(qiáng)迫函數(shù)為 MN 12P , a ecm2m2iim i m1i 1mm(3-3) d MN12Q , a m m1 e cm m2b2mmeiii i i 1mm在變換空間中，這些函數(shù)寫成為 Qx (3-4)Py Qy y 2y y J2 yy ，而這里 J 是變換的 Jacobi 矩陣行列式的值，J = x y , x x y y , x 2 y 222 和代數(shù)法相比，偏微分方程法的計(jì)算較為復(fù)雜，需要求解偏微分方程，右端強(qiáng)迫函數(shù)的選取也不方便，可以通過(guò)強(qiáng)迫函數(shù)的選取得到希望的網(wǎng)格。213.2.3 超單元法等參超單元網(wǎng)格生成

9、法是由 Zienkiewicz 和31,42，其基本Phillips 最早是將劃分區(qū)域分為自然坐標(biāo)系19成更簡(jiǎn)單的子塊（超單元），子塊被5中的直角形（正方形），根據(jù)每個(gè)方向給定的分級(jí)權(quán)系數(shù)，將自然坐標(biāo)系中的直角形（正方形）離散并反向變換回超單元。在三維情況使用超單元生成法時(shí)，要?jiǎng)澐值膶?shí)體用 20 節(jié)點(diǎn)六面體單元分成數(shù)個(gè)子區(qū)域（圖3-3）。在超單元內(nèi)一點(diǎn)的坐標(biāo) x，y 和 z 和自然坐標(biāo)之間的關(guān)系由下式確定163202020 x Ni xi , y Ni yi , z Ni zi1(3-5)92i 1i1i 1式中，Ni(, , )是與按曲線坐標(biāo)系， 和定義的與每個(gè)節(jié)點(diǎn)相聯(lián)系的形函數(shù)， 和的取值

10、范圍為1 到 1。 20 節(jié)點(diǎn)六面體單元形函數(shù)在各角點(diǎn)定義為圖3-3 20 節(jié)點(diǎn)六面體單元Fig. 3-3 20-node hexahedral element, i 1, 2, , 8(3-6)在各棱邊的中間節(jié)點(diǎn)定義為N ( , , ) 1 1 2 1 1 , i 9,11, 17, 19iii41 1N ( , , ) 1 1 2, i 10,12, 18, 20(3-7)iii4N ( , , ) 1 1 1 1 2 , i 13,14, 15,16iii4如節(jié)點(diǎn)的坐標(biāo)已知，則與特定坐標(biāo)， 和相對(duì)應(yīng)的點(diǎn)的坐標(biāo)可由式(3-5)計(jì)算出。沿每個(gè)局部自然坐標(biāo)，給定劃分的數(shù)目和劃分的間距，這些子區(qū)

11、域（超單元）被劃分為更小的單元。以為例，劃分時(shí)自然坐標(biāo)按下式增加 i (3-8)0 = 1。法計(jì)算簡(jiǎn)單，通過(guò)對(duì)自然坐標(biāo)的不同劃分，可以得到密度變化的網(wǎng)格，但超單元這里W超單元描述的邊界和劃分區(qū)域的實(shí)際邊界存在誤差，不能適應(yīng)復(fù)雜的邊界。223.3 非結(jié)構(gòu)網(wǎng)格四邊形/六面體網(wǎng)格生成法由于算法的特點(diǎn)，自動(dòng)非結(jié)構(gòu)網(wǎng)格生成算法主要是三角形或四面體網(wǎng)格劃分，有許多文獻(xiàn)討論了三角形或四面體網(wǎng)格劃分方法和相關(guān)的技術(shù)43,44,45,46,47,48,49,50,51,52,53,54，許多也提供了三角形或四面體網(wǎng)格劃分功能55。非結(jié)構(gòu)網(wǎng)格生成的特點(diǎn)是能適應(yīng)復(fù)雜的邊界，需要的人工干預(yù)少，計(jì)算費(fèi)時(shí)。Steven.

12、 J. Owen55總結(jié)了當(dāng)前學(xué)術(shù)界和工業(yè)界使用的網(wǎng)格劃分方法和可以利用的。非結(jié)構(gòu)六面體網(wǎng)格生成方法主要有法、轉(zhuǎn)換法、鋪層算法和須段編織算法等。不象三角形網(wǎng)格劃分算法向四面體網(wǎng)格劃分算法擴(kuò)展那樣直接，從二維四邊形網(wǎng)格劃分算法向三維六面體網(wǎng)格算法的擴(kuò)展并不是直接的。3.3.1法當(dāng)劃分區(qū)域的幾何特性合適時(shí)，四邊形或六面體網(wǎng)格生成法可以得到相當(dāng)好的網(wǎng)格56。雖然法被認(rèn)為是結(jié)構(gòu)化網(wǎng)格生成法，但許多非結(jié)構(gòu)網(wǎng)格劃分都提供這種方法。當(dāng)應(yīng)用映射法時(shí)，要求劃分區(qū)域的對(duì)邊具有相同的子劃分。在三維情況下，將要?jiǎng)澐謪^(qū)域的拓?fù)浣Y(jié)構(gòu)的立方體的每個(gè)相對(duì)的面必須有相同的表面網(wǎng)格。這種情況在幾何形狀復(fù)雜是不可能達(dá)到，或者需要人

13、工干預(yù)將幾何區(qū)域劃分成可的子區(qū)域并分配適當(dāng)?shù)拈g隔。等57,58使用的方法是將劃分區(qū)域分成子區(qū)域，每個(gè)子區(qū)域?qū)?yīng)于不同的模塊，每個(gè)模塊中的網(wǎng)格已預(yù)先劃分，模塊中的網(wǎng)格為子區(qū)域中的實(shí)際網(wǎng)格。這里模塊中的網(wǎng)格可以是結(jié)構(gòu)化的，也可以是非結(jié)構(gòu)化的，但要求相鄰子區(qū)域的交接面上的表面網(wǎng)格具有相同的拓?fù)浣Y(jié)構(gòu)。該方法的特點(diǎn)是可以針對(duì)子區(qū)域的幾何特性和物理特性選擇相應(yīng)的模塊，生成的網(wǎng)格質(zhì)量有 ANSYS（ANSYS Inc.）、CUBIT（Sandia高，缺點(diǎn)是通用性差。提供National Laboratories）、goe Inc.）等。法生成網(wǎng)格的Wrap（Raindrop GeoInc.）、Hyperme

14、sh（Altair Computing掃描線目標(biāo)源圖3-4 掃描法生成六面體單元Fig. 3-4 Hex elements generated by s）網(wǎng)格生成法是另一類六面體網(wǎng)格生成方法59,60,61，有時(shí)被認(rèn)為是二掃描法（S維半網(wǎng)格生成方法。該方法是將四邊形網(wǎng)格沿空間一條曲線掃描，給定一定的間隔，使用該四邊形網(wǎng)格的拓?fù)浣Y(jié)構(gòu)生成給定層的六面體網(wǎng)格。這種方法歸納為給定源表面和目標(biāo)表面劃分某些類型的體積區(qū)域的技術(shù)。假定源表面和目標(biāo)表面具有相同的拓?fù)浣Y(jié)構(gòu)并被可劃分表面相連接，23源表面上的四邊形網(wǎng)格掃過(guò)體積區(qū)域生成六面體網(wǎng)格（圖3-4）。提供掃描法生成網(wǎng)格的有ANSYS（ANSYS Inc.）

15、、CUBIT（Sandia National Laboratories）、Hypermesh（Altair Computing Inc.）等。3.3.2 非結(jié)構(gòu)四邊形網(wǎng)格劃分非結(jié)構(gòu)四邊形網(wǎng)格劃分算法一般可分為兩類：直接方法和間接方法。間接方法是先將劃分區(qū)域劃分為三角形，然后采用各種算法將三角形轉(zhuǎn)換為四邊形。直接法是在劃分區(qū)域上直接生成四邊形，不需要先生成三角形。圖3-5 三角形分為四邊形生成的網(wǎng)格Fig. 3-5 Quad mesh generated by splitting each triangleo three quads圖3-6 三角形組合為四邊形生成的網(wǎng)格Fig. 3-6 Quad

16、-dominant mesh generated by combining triangles1三角形分解/合并成四邊形（間接法）最簡(jiǎn)單的間接四邊形網(wǎng)格生成是將每個(gè)三角形劃分為三個(gè)四邊形，如圖3-5所示，這種方法保證可以生成全四邊形的網(wǎng)格，但網(wǎng)格中大量存在的不規(guī)則節(jié)點(diǎn)導(dǎo)致低的網(wǎng)格質(zhì)量。另一種算法是組合相鄰的兩個(gè)三角形形成一個(gè)四邊形，如圖3-6所示，這種方法可以增加網(wǎng)格的質(zhì)量，但同時(shí)會(huì)o62給出的算法留下大量的三角形單元。在三角形組合次序上采取某些措施，可以改進(jìn)該方是通過(guò)分析可能的合并，按照合并后四邊形的優(yōu)劣結(jié)果決定合并的次序，得到一個(gè)包含最少三角等63對(duì)Lo 的算法進(jìn)行了改進(jìn)。Johnson

17、等64提出采用增加形單元的四邊形占優(yōu)的網(wǎng)格。閔對(duì)局部單元的劃分和交換的辦法來(lái)提高四邊形單元的數(shù)量和質(zhì)量。Lee 等65通過(guò)引入局部三角形的分解擴(kuò)展了 Lo 的算法62，并使用了波前的概念，如果區(qū)域邊界上初始的單元邊數(shù)是偶數(shù)，這種方法可以組全四邊形網(wǎng)格。Owen 等66Q-Morph 算法用于從三維的三角形網(wǎng)格中生24成四邊形網(wǎng)格，該算法采用決定三角形轉(zhuǎn)換的序列，使用三角形網(wǎng)格中已存在的邊或額外的節(jié)點(diǎn)或?qū)植咳切芜M(jìn)行變換來(lái)形成四邊形，任意數(shù)量的三角形可以合并生成一個(gè)四邊形。Q-Morph 算法可以生成沿區(qū)域邊界排列的很好的四邊形，但在區(qū)域還是存在數(shù)目有限的不規(guī)則節(jié)點(diǎn)。Q-Morph 算法在有限

18、元ANSYS 上得到了實(shí)現(xiàn)。由于三角形分解合并為四邊形的所有操作都限于局部，故間接算法的速度非?？欤恍枰笾苯铀惴菢舆M(jìn)行相交檢查。其缺點(diǎn)是網(wǎng)格中會(huì)留下許多不規(guī)則的節(jié)點(diǎn)。對(duì)某些應(yīng)用而言，即使留下的不規(guī)則節(jié)點(diǎn)是少數(shù)，也不能保證沿區(qū)域邊界的單元能以希望的形式排列。可以通過(guò)拓?fù)浣Y(jié)構(gòu)的調(diào)整去掉這些不規(guī)則的節(jié)點(diǎn)從而提高網(wǎng)格的質(zhì)量。2柵格法Baehmann 等67采用的四叉樹分解方法是將目標(biāo)區(qū)域用盡量小的正方形圈定，然后把這個(gè)正方形分解為四個(gè)子區(qū)域，對(duì)每個(gè)子區(qū)域測(cè)試其是否完全落在目標(biāo)區(qū)域之外，或是否滿足密度控制要求，若滿足則停止對(duì)區(qū)域的細(xì)分，否則細(xì)分下去，遞歸執(zhí)行直到到達(dá)預(yù)定的要求，留下合適的四邊形單元

19、，調(diào)節(jié)節(jié)點(diǎn)位置保持與邊界的一致。柵格法（Grid Method）68,69是將一包含目標(biāo)區(qū)域的柵格放置在目標(biāo)區(qū)域上面，除去落在區(qū)域外的柵格單元，并將與區(qū)域邊界相交的柵格單元進(jìn)行調(diào)整或裁剪，然后在柵格邊界單元和區(qū)域邊界之間生成四邊形單元（圖3-7）。區(qū)域的網(wǎng)格最終的網(wǎng)格圖3-7 基于柵格的四邊形網(wǎng)格生成Fig. 3-7 Quadrilateral mesh generated by rigid-based methodTalbert 等70使用的分解方法是將區(qū)域迭代劃分為簡(jiǎn)單的多邊形子區(qū)域，這些子多邊形的形狀與有限數(shù)量的模板相匹配，然后在這些多邊形中四邊形單元。Tam 和 Armstrong71

20、利用區(qū)域的中間軸（Medial Axis）將區(qū)域分解為簡(jiǎn)單凸多邊形區(qū)域，然后使用一套模板將四邊形網(wǎng)格區(qū)域。平面對(duì)象的中間軸是對(duì)象中可放入最大圓的圓心軌跡，即沿對(duì)象滾動(dòng)的最大直徑的圓的圓心軌跡（圖3-8）。中間軸由數(shù)個(gè)直線或曲線段組成，采用一定的規(guī)則用這些線段將區(qū)域分解為凸的子區(qū)域，再采用界拓?fù)湎嗳莸木W(wǎng)格。的方法在子區(qū)域上生成與子區(qū)域邊有限四叉樹分解法和基于網(wǎng)柵的分解法在區(qū)域可以生成質(zhì)量很高的單元，但在區(qū)域邊界上要滿足邊界一致性要求，邊界單元的質(zhì)量較差，不適合塑性成形變形大、網(wǎng)格容易發(fā)生畸變的應(yīng)用領(lǐng)域。分解為子區(qū)域的劃分方法，通過(guò)控制所劃分的子區(qū)域的形狀和增加可模板的數(shù)量，25可以提高子區(qū)域和模

21、板的相似性，從而使網(wǎng)格質(zhì)量得到提高，然而對(duì)于一般幾何形狀的區(qū)域，難于做到通用性。中間軸圖3-8 二維對(duì)象的中間軸Fig. 3-8 Medial axis of a 2D object圖3-9 鋪路算法示意Fig. 3-9 Demo of paving algorithm3四邊形網(wǎng)格生成的Zhu 等72最早提出了使用的四邊形網(wǎng)格的生成算法。該算法是先在邊界上生成初始節(jié)點(diǎn)，然后向區(qū)域投射邊段，先使用傳統(tǒng)的三角形網(wǎng)格產(chǎn)生兩個(gè)三角形單元，再將兩個(gè)三角形組合為一個(gè)四邊形。Blacker 和 Stephenson鋪路算法73（Paving Algorithm）是重復(fù)地在區(qū)域邊界放一等74將鋪路層或者說(shuō)鋪一

22、排單元，算法從區(qū)域邊界開始向形成四邊形單元（圖3-9）。Cass算法推廣到了三維表面的劃分，White 和 Kinney 重新設(shè)計(jì)了鋪路算法75，單元不是生成一排，而是一個(gè)接一個(gè)的生成。方法。CUBIT（Sandia National Laboratories）等提供了鋪路算法的網(wǎng)格劃分鋪路算法能生成邊界上具有較高質(zhì)量的網(wǎng)格，是一種較好的非結(jié)構(gòu)網(wǎng)格生成方法。3.3.3 非結(jié)構(gòu)六面體網(wǎng)格生成類似于非結(jié)構(gòu)四邊形網(wǎng)格生成方法，非結(jié)構(gòu)六面體網(wǎng)格生成也可以分為直接法和間接法。圖3-10 四面體分解為六面體ition of a tetrahedrono four hexahedral elementsFi

23、g. 3-10261四面體分解為六面體76（間接法）如果體積區(qū)域可以劃分為四面體網(wǎng)格，則每個(gè)四面體單元可以劃分為四個(gè)六面體單元，如圖77采用 103-10所示。一般情況下這種方法生成的單元質(zhì)量較差，在有限元分析時(shí)很少采用。節(jié)點(diǎn)的四面體單元分解為四個(gè)六面體單元的方法，可以提高邊界處單元與劃分區(qū)域邊界的一致性，但對(duì)單元質(zhì)量并沒(méi)有實(shí)質(zhì)性的提高。還沒(méi)有文獻(xiàn)對(duì)四面體單元組。2基于柵格的方法六面體單元的方法進(jìn)行基于柵格的六面體網(wǎng)格生成方法包括（有限）八叉樹分解法78和正則網(wǎng)格覆蓋法79。這里先采用八叉樹分解或網(wǎng)格覆蓋的方法，在劃分區(qū)域生成合適的三維六面體單元網(wǎng)格，然后在邊界處加入或調(diào)整六面體網(wǎng)格，將規(guī)則網(wǎng)

24、格與邊界之間的空白區(qū)域充滿，并滿足邊界一致性要求。基于柵格的方法在所劃分的體積區(qū)域的邊界處生成的單元質(zhì)量不高，通常單元沿邊界的排列也不好，同時(shí)生成的網(wǎng)格取決于成的單元除邊界處外大小相同，具有明顯單元尺寸過(guò)渡的網(wǎng)格。于柵格的六面體網(wǎng)格生成方法。六面體網(wǎng)格的取向，即與坐標(biāo)軸的方向有關(guān)。正則網(wǎng)格法生等80和Schneider81對(duì)八叉樹算法進(jìn)行了修改，可以生成Men （MARCysis Research Corpoion）等提供了基3中間面法（Medial Surface Method）中間面法82,83,84是劃分四邊形網(wǎng)格中間軸法在三維的直接擴(kuò)展。這里劃分區(qū)域被分解為由中間面分開的子區(qū)域的集合，

25、中間面可以認(rèn)為是在區(qū)域滾過(guò)的最大球的球心點(diǎn)組成的面。由中間面分開所形成的子區(qū)域形成可劃分網(wǎng)格的區(qū)域，這些區(qū)域?qū)?yīng)于一系列的具有所期望的拓?fù)浣Y(jié)構(gòu)的模板，整個(gè)區(qū)域用后的六面體網(wǎng)格充滿，使用線性編程來(lái)保證從一個(gè)區(qū)域到另一個(gè)區(qū)域單元?jiǎng)澐种g的匹配。這種方法適合于某些具有特定幾何形狀的區(qū)域的劃分，但對(duì)一般幾何形狀的區(qū)域可靠性差，生成中間面及規(guī)定中間面定義的所有區(qū)域還是一個(gè)難于解決。中間面劃分算法在 CADfix（FEGS.）、TurboMesh（SolidPo）等中得到實(shí)現(xiàn)。4鋪層算法（Plastering Algorithm）鋪層算法85是鋪路算法在三的拓展，這種算法是先將單元放在邊界處，然后逐層向劃

26、分區(qū)域中心推進(jìn)（圖3-11），定義一套探索程序決定單元形成的次序。和其它推進(jìn)類似，所有的四邊形當(dāng)前的波前，每個(gè)的四邊形向體積區(qū)域投射形成六面體。此外，算法必須判斷與表面的相交、何時(shí)和如何連接原來(lái)存在的節(jié)點(diǎn)來(lái)縫合留下的空間。隨著算法的繼續(xù)，在區(qū)域可能留下復(fù)雜的空隙，在某些情況下這些區(qū)域不能用全六面體單元填充。為了方便地向內(nèi)放入六面體，鋪層算法有時(shí)必須修改已生成的單元。鋪層算法在大規(guī)模問(wèn)題上應(yīng)用的可靠性還沒(méi)有被證明。雖然可以成功地在邊界上放上數(shù)層六面體單圖3-11 鋪層法生成單元的過(guò)程Fig. 3-11 Plastering Pros forming elements at the boundar

27、y元，但相交和封閉算法的可靠性不高。等提供了鋪層算法劃分網(wǎng)格的方法。CUBIT275須段編織算法Tautges 等86須段編織算法（Whisker Weaving Algorithm）是基于空間纏連集（SpatialTwist Continuum，STC）的一種拓?fù)鋭澐址椒āTC 是建立在全六面體網(wǎng)格對(duì)等體上的一個(gè)整體構(gòu)造集，在 STC 中，用相交面的排列來(lái)對(duì)等地表示網(wǎng)格的實(shí)體，這些面在各個(gè)方向上平分六面體單元并形成相交的曲線和頂點(diǎn)，STC 封裝了全六面體網(wǎng)格的約束信息。平分每個(gè)單元的面有三個(gè)，這些面稱為須段面（Whisker Sheet），可以認(rèn)為它是一個(gè)一般的的表面，這個(gè)表面在網(wǎng)格中延續(xù)

28、，在相對(duì)的兩個(gè)面的中間平分每個(gè)六面體。圖3-12是由兩個(gè)單元組成的網(wǎng)格，它有 4 個(gè)須段面，編號(hào) 1 至 4，這些面的排列就稱為 STC。一個(gè)須段面代表一層相鄰的單元之間具有共享面的單元。兩個(gè)須段面相交形成的曲線稱為須段弦（Whisker Chord），圖中有 5 個(gè)須段弦，用 A 至 E 標(biāo)識(shí)，須段弦表示網(wǎng)格中相鄰單元之間具有共享面的一列單元，如須段線 E 定義了這個(gè)網(wǎng)格中兩個(gè)單元組成的單元列。每個(gè)六面體單元對(duì)應(yīng)于一個(gè)對(duì)等的稱為 STC 頂點(diǎn)（STC Vertex）的實(shí)體，它是三個(gè)須段面相交形成的交點(diǎn)（同樣是三條須段線的交點(diǎn)）。須段面與幾何表面的交線稱為環(huán)（Loop），它是一條封閉曲線，表示

29、一圈表面四邊形，每個(gè)相鄰四邊形具有共享的邊。BD12aAbC34圖3-12 兩個(gè)單元的 STCFig. 3-12 The STC for a solid comed of two hexahedra須段面是空間中和相交的二維表面，它們之間的相交的幾何位置是隱含的，在須段編織中，須段面用平面面圖（Sheet Diagram）表示（圖3-13）。當(dāng)須段面變平時(shí)，它的環(huán)形成一個(gè)多邊形，它的須段弦是多邊形內(nèi)的線段。須段面上的須段弦在面圖上用一個(gè)邊列表示，為與須段弦相區(qū)別稱為面弦（Sheetchord）。須段編織算法從定義一個(gè)封閉三維體積實(shí)體的四邊形表面網(wǎng)格開始的，因此必須先生成劃分區(qū)域表面的四邊形網(wǎng)格

30、。須段編織算法分為四步構(gòu)造初始的環(huán)、須段面、須段弦和面弦。從給定的表面網(wǎng)格取一未遍歷過(guò)的邊，從這條邊開始，向四邊形的對(duì)邊遍歷，直到從另一邊回到開始邊就形成了一個(gè)環(huán)。當(dāng)所有的四邊形都被兩個(gè)環(huán)穿過(guò)時(shí)所有的環(huán)都形成了。環(huán)生成后，構(gòu)造每個(gè)須段面的須段弦和對(duì)應(yīng)的面弦。找出三條兩兩相鄰須段弦。28E在須段編織算法中，三條須段弦相交形成 STC 頂點(diǎn)組成六面體，在面圖上表現(xiàn)為三張須段面上三對(duì)相交的面弦。使三條須段弦相交形成一個(gè)新的須段六面體。檢查連接的有效性，連接須段弦。重復(fù)步驟（2）（4）直到完成編織。完成的標(biāo)志是沒(méi)有懸掛的須段留下。完成的編織表示全六面體網(wǎng)格的整體的連接性，而須段編織并沒(méi)有產(chǎn)生六面體單元

31、幾何位置的任何信息，環(huán)須段弦也沒(méi)有決定體積區(qū)域?qū)嶋H單元或節(jié)點(diǎn)的位置，這些位置信息可以通過(guò)使用初始構(gòu)造算法（Primal Construction Algorithm）得到。須段編織算法可以生成質(zhì)量較高的網(wǎng)格，但對(duì)范圍廣泛的問(wèn)題算法的穩(wěn)定性（Robustness）和可靠性還沒(méi)有得到證實(shí)。6六面體單元占優(yōu)網(wǎng)格圖3-13 須段面圖Fig. 3-13 A whisker sheet diagram由于大多數(shù)六面體網(wǎng)格生成方法的穩(wěn)定性較差，有些研究者建議使用六面體四面體混合網(wǎng)格。在許多情況下，六面體占優(yōu)的網(wǎng)格劃分的方法是令人滿意的。Owen 給出的一個(gè)簡(jiǎn)單方法是手動(dòng)劃分區(qū)域成可劃分的子區(qū)域和幾何上較復(fù)雜

32、的子區(qū)域，在復(fù)雜的子區(qū)域上定義四面體單元，在四面體單元和六面體單元的的交界處使用過(guò)渡單元。在 ANSYS 網(wǎng)格生成這種算法劃分網(wǎng)格。中可以選擇Tuchinsky 等87,88使用的方法是組合鋪層算法和四面體網(wǎng)格劃分算法，使用鋪層算法盡可能地將六面體單元推進(jìn)入劃分體積區(qū)域，在區(qū)域內(nèi)留下的空隙用四面體單元填充，也可選擇五面體單元作為六面體單元和四面體單元的過(guò)渡。 CUBIT 可以選擇這種方法生成六面體占優(yōu)的網(wǎng)格。3.4 表面網(wǎng)格生成當(dāng)前許多網(wǎng)格生成與任意表面單元的生成有關(guān)。典型地，這些表面是商業(yè) CAD包生成的 NURBS 曲面。表面單元可以直接作為殼單元使用，也可作為體積網(wǎng)格的輸入數(shù)據(jù)。將二維網(wǎng)

33、格生成算法推廣到三維表面網(wǎng)格的劃分需要作某些修改。表面網(wǎng)格生成算法可分為參數(shù)空間法和直接劃分法兩類。3.4.1 參數(shù)空間法所有的 NURBS 曲面都是基于二維參數(shù)空間(u, v)表示的，而二中的網(wǎng)格劃分通常效率較高。參數(shù)空間算法（Parametric Space Algorithm）是先在劃分表面的二維參數(shù)空間(u, v)中生成網(wǎng)格，然后將參數(shù)空間的坐標(biāo)(u, v)回實(shí)際坐標(biāo)(x, y, z)。這種方法的缺點(diǎn)是在參數(shù)空間中生成的單元回劃分表面時(shí)，并不總是形成三中的好形狀的單元?？梢圆扇煞N措施來(lái)解決這個(gè)問(wèn)題：修改或重新參數(shù)化曲面的參數(shù)表示，以得到從參數(shù)空間到實(shí)際空間的合適的；修改網(wǎng)格生成算法，

34、使在二維區(qū)域劃分的扁的或不等向的單元能元。為三維好形狀、等向的單29為了得到一個(gè)好的參數(shù)化表示，要求在定義域內(nèi)表面偏導(dǎo)數(shù)和(u, v)變化不是太大。文獻(xiàn) 89中有精確弧長(zhǎng)的重新參數(shù)化定義，但花費(fèi)過(guò)大。通過(guò)有選擇計(jì)算區(qū)域內(nèi)的表面偏導(dǎo)數(shù)，調(diào)節(jié)局部的(u, v)的值保持(u, v)的大小大致為常數(shù)，可以定義近似的弧長(zhǎng)參數(shù)化或稱翹曲參數(shù)空間（Warped Parametric Space）。在許多情況下，翹曲參數(shù)空間能生成合理的網(wǎng)格，但還有許多問(wèn)題是重新參數(shù)化不能適當(dāng)解決的，由于這種原因，多數(shù)有關(guān)表面網(wǎng)格劃分的文獻(xiàn)集中了第二種方法：在二中形成不等向的網(wǎng)格，然后成三的等向網(wǎng)格。實(shí)踐中使用的一個(gè)通用的方法

35、是利用表面偏導(dǎo)數(shù)、u 和v 這些可以方便地從 NURBS 曲面中計(jì)算出的量。Gee 和 Borouchaki55使用了一個(gè)由表面第一基本形式導(dǎo)出的表面度量，該度量是一個(gè) 22 矩陣，用于變換參數(shù)空間的向量和距離。Cuillire90使用節(jié)點(diǎn)密度優(yōu)先的波前算法完成了參數(shù)表面的三角化網(wǎng)格劃分，其波前是在劃分表面的邊界上形成，然后反回參數(shù)空間，在參數(shù)空間中使用通常的波前算法將網(wǎng)格劃分，劃分時(shí)同樣采用了表面度量來(lái)決定新點(diǎn)的位置。Lee 和 Joun 給出了一種基于的四邊形表面網(wǎng)格生成方法91，這種算法是將復(fù)雜的、不利于直接的表面裁剪成三角形，三角形邊界上的一個(gè)點(diǎn)為參數(shù)區(qū)域中的一條線段，為避免出現(xiàn)此種不

36、利情形，將區(qū)域劃分四邊形網(wǎng)格，然后反MARC/Men、ANSYS 等參數(shù)區(qū)域再次，成為三角形的再區(qū)域，在再回初始區(qū)域的實(shí)際網(wǎng)格。中提供了參數(shù)空間法生成網(wǎng)格的算法。3.4.2 直接劃分法直接法三維表面網(wǎng)格生成方法是直接在劃分的幾何區(qū)域上形成網(wǎng)格而不管幾何區(qū)域的參數(shù)表示。在沒(méi)有參數(shù)化表示可用或表面參數(shù)化非常差時(shí)，可以使用直接法。Lau 和 Lo92使用波前算法在任意空間曲面生成了三角形網(wǎng)格。他們的方法是給定一個(gè)波前的線段，檢查與邊界上或區(qū)域內(nèi)的的一個(gè)節(jié)點(diǎn)生成新三角形的可能性，根據(jù)曲面曲率、單元質(zhì)量和給定的單元密度分布函數(shù)來(lái)優(yōu)化生成的單元，為保持新生成的節(jié)點(diǎn)位于曲面之上，需要將生成的節(jié)點(diǎn)投射到曲面上

37、。Cass 等74通過(guò)引入節(jié)點(diǎn)投射、決定波前網(wǎng)格 、支持周期性表面和單元的平滑等技術(shù)，將鋪路算法擴(kuò)展到了三維曲面的四邊形網(wǎng)格生成。在CUBIT 中可以利用該方法生成表面四邊形網(wǎng)格。3.5 網(wǎng)格后處理一般來(lái)說(shuō)，無(wú)論用何種方法生成網(wǎng)格，總會(huì)有一些形狀不好的單元產(chǎn)生，如內(nèi)角太大的單元、長(zhǎng)條單元、薄單元等。這些單元的存在本身會(huì)引起有限元解答的精度降低，甚至?xí)?dǎo)致解的不收斂。在金屬塑性成形有限元模擬過(guò)程中，這些形狀不好的單元會(huì)很快發(fā)生，導(dǎo)致頻繁的網(wǎng)格再劃分。為了減少模擬時(shí)間，提高模擬精度，必須采用某種或幾種方法，將這些性狀不好的單元改性以改進(jìn)網(wǎng)格的整體質(zhì)量。對(duì)六面體單元，可用于判斷單元質(zhì)量的數(shù)據(jù)有單元的

38、長(zhǎng)寬比，即單元的最大與最小邊長(zhǎng)的比值。每個(gè)節(jié)點(diǎn)球面度的誤差，球面度的定義參見(jiàn) 6.4 節(jié)?？梢匀∶總€(gè)節(jié)點(diǎn)的球面度與標(biāo)準(zhǔn)值差值的百分比絕對(duì)值作為判斷標(biāo)準(zhǔn)。體積測(cè)量。六面體每個(gè)節(jié)點(diǎn)連接三條邊，以三條邊為起始于該節(jié)點(diǎn)三個(gè)向量 r1，r2 和 r3，30對(duì)每個(gè)節(jié)點(diǎn)計(jì)算 V = (r1r2)r3，采用 C = |V Ve| / Ve 作為體積測(cè)量的判斷標(biāo)準(zhǔn)，Ve 為單元的體積。4單元面的翹曲度。六面體單元每個(gè)面由四個(gè)節(jié)點(diǎn)組成，有四條邊，每三個(gè)節(jié)點(diǎn)確定一個(gè)平面，平面與另一個(gè)節(jié)點(diǎn)的距離為 d，單元面的翹曲度定義為最大的 d 與最短邊長(zhǎng)的比值。5單元的度。六面體單元由三對(duì)相對(duì)的面，單元的度定義為相對(duì)面夾角的銳

39、角。改進(jìn)網(wǎng)格質(zhì)量的方法主要有兩類：節(jié)點(diǎn)平滑和網(wǎng)格整理。節(jié)點(diǎn)平滑是在不改變網(wǎng)格拓?fù)浣Y(jié)構(gòu)的條件下，通過(guò)調(diào)整節(jié)點(diǎn)的位置來(lái)提高網(wǎng)格質(zhì)量。網(wǎng)格整理通常是指通過(guò)合并單元、劃分單元、單元節(jié)點(diǎn)的重新分派等改變網(wǎng)格拓?fù)浣Y(jié)構(gòu)的方法改善網(wǎng)格質(zhì)量的方法，由于網(wǎng)格結(jié)構(gòu)的特點(diǎn)，六面體網(wǎng)格很少采用網(wǎng)格整理的方法來(lái)提高網(wǎng)格的質(zhì)量。節(jié)點(diǎn)平滑處理技術(shù)大多數(shù)節(jié)點(diǎn)平滑程序是采用迭代過(guò)程將節(jié)點(diǎn)重定位來(lái)改善局部單元質(zhì)量的。有多種平滑技術(shù)已被提出，通常這些方法可分為平均方法、基于優(yōu)化的方法、基于物理的方法、中間節(jié)點(diǎn)定位等四種方法。平均方法在所有平滑算法中，Laplace 平滑算法是應(yīng)用最廣的一種，這種算法是把網(wǎng)格的一個(gè)節(jié)點(diǎn)放在與該節(jié)點(diǎn)用邊

40、相連的節(jié)點(diǎn)的平均位置上。如圖3-14a 所示，每個(gè)節(jié)點(diǎn)的新位置由與該節(jié)點(diǎn)相連節(jié)點(diǎn)的平均坐標(biāo)決定1nP P(3-9)ini1這里 P 是要平滑的節(jié)點(diǎn)的位置向量，Pi 直接連接到節(jié)點(diǎn) P 的鄰近的節(jié)點(diǎn)的位置向量，n 為連接到該節(jié)點(diǎn)的節(jié)點(diǎn)數(shù)。Laplace 平滑是共享邊的平滑，而等參平滑(Isoparametric Smoothing)是共享面的平滑，如圖 3-14b，使用所有與節(jié)點(diǎn) PS 相連的面，直接與 PS 相連的節(jié)點(diǎn)坐標(biāo)取為正，而在面上對(duì)角的節(jié)點(diǎn)坐標(biāo)取為減。圖中括號(hào)內(nèi)的加減號(hào)數(shù)目表示同一節(jié)點(diǎn)的加減次數(shù)。經(jīng)過(guò)小量的修改，這些方法可應(yīng)用于任何形狀的單元。絕大多數(shù)這樣的平滑程序需要在網(wǎng)格中所有節(jié)點(diǎn)

41、進(jìn)行多次迭代，直到每個(gè)節(jié)點(diǎn)移動(dòng)的距離小于給定的偏差。雖然這些方法存在某些問(wèn)題，但程序?qū)崿F(xiàn)簡(jiǎn)單，因而得到了廣泛應(yīng)用。類似于 Laplace 節(jié)點(diǎn)平滑方法，基于周圍環(huán)繞。Canann 等93對(duì)常用節(jié)點(diǎn)優(yōu)節(jié)點(diǎn)和單元幾何特性的化方法進(jìn)行了總結(jié)。平均，產(chǎn)生了多種迭代節(jié)點(diǎn)重單獨(dú)使用平均方式平滑節(jié)點(diǎn)有時(shí)會(huì)降低局部單元的質(zhì)量，為避免降低單元質(zhì)量的現(xiàn)象出現(xiàn)，平均方法通常也要對(duì)節(jié)點(diǎn)的移動(dòng)施加某種形式的約束，對(duì)節(jié)點(diǎn)移動(dòng)前后局部單元的質(zhì)量進(jìn)行比較，只有單元質(zhì)量有所改善時(shí)才移動(dòng)節(jié)點(diǎn)的位置。這些方法通常稱為約束Laplace 平滑（ConstrainedLaplan Smoothing）。31P3()P4()P6P2()

42、P3P4PP5()SP1()P13()PP10()P6()P ()P9()8PP1P52P7()P ()P12()11（b）（a）圖3-14 Laplace 平滑和等參平滑F(xiàn)ig. 3-14 Laplan smoothing and isoparametric smoothing基于優(yōu)化的方法基于優(yōu)化的平滑技術(shù)測(cè)量環(huán)繞一個(gè)節(jié)點(diǎn)單元的質(zhì)量，通過(guò)計(jì)算單元質(zhì)量關(guān)于節(jié)點(diǎn)位置的局部梯度來(lái)優(yōu)化單元的質(zhì)量，節(jié)點(diǎn)沿梯度增加的方向移動(dòng)，直到達(dá)到優(yōu)化的目標(biāo)94,77。單元質(zhì)量好壞由單元長(zhǎng)度方向的相對(duì)關(guān)系和單元各個(gè)面（邊）之間的相對(duì)關(guān)系來(lái)衡量，反映在單元變形上與拉壓變形和剪切變形有關(guān)，而單元的 Jacobi 矩陣中

43、的各項(xiàng)反映了這些變形，故常用 Jocobi 矩陣分量來(lái)衡量單元的質(zhì)量。有關(guān)的算法如下將單元的 Jacobi 矩陣正則化，去掉單元尺寸的影響J ijJ (3-10)ij1Jn這里 n 是問(wèn)題的維數(shù)。2引入能描述單元畸變和應(yīng)變之間模擬關(guān)系的 Green 應(yīng)變，它的一項(xiàng)可表示為 Jacobi 矩陣的函數(shù)Cij Jki J kj(3-11)3單元畸變的度量D C C 1 C2(3-12)ij ijkkn它是 Jacobi 矩陣的函數(shù)，也即是節(jié)點(diǎn)坐標(biāo)的函數(shù)。對(duì)整個(gè)區(qū)域得到一個(gè)畸變的度量值，以此作為目標(biāo)函數(shù)進(jìn)行優(yōu)化。目標(biāo)函數(shù)變量較多，可采用共軛梯度法(Conjugate Gradient Method)等

44、進(jìn)行優(yōu)化，得到優(yōu)化的節(jié)點(diǎn)坐標(biāo)值。除此之外，還有多種關(guān)于單元質(zhì)量度量的分析表達(dá)可用于單元質(zhì)量的優(yōu)化，這些量包括單元32每條邊連接的兩個(gè)面的夾角、每個(gè)節(jié)點(diǎn)處的球面度（對(duì)二維單元為平面角）、單元的扁率（ElementAspect Ratio ）、單元面積（周長(zhǎng)）與體積（面積）的比值等。所有這些分析量可以向上面所述那樣單獨(dú)使用，但由于單元和網(wǎng)格的質(zhì)量不是任意一個(gè)分析量可以完全描述的，故一般的做法是將幾個(gè)分析量采用組合的方法，組成一個(gè)目標(biāo)函數(shù)，對(duì)這個(gè)目標(biāo)函數(shù)進(jìn)行優(yōu)化。一般來(lái)說(shuō)目標(biāo)函數(shù)的自變量是網(wǎng)格節(jié)點(diǎn)的坐標(biāo)，變量較多，可采用適合多變量?jī)?yōu)化的方法進(jìn)行優(yōu)化，得到優(yōu)化的節(jié)點(diǎn)坐標(biāo)?；趦?yōu)化的方法可以得到質(zhì)量?jī)?yōu)良

45、的網(wǎng)格，但由于變量多、約束條件多，導(dǎo)致計(jì)算量過(guò)大且不易收斂。有的研究者95建議將 Laplace 節(jié)點(diǎn)平滑技術(shù)和基于優(yōu)化的方法組合使用來(lái)優(yōu)化網(wǎng)格，通常大部分時(shí)間采用 Laplace 平滑，僅3.5.1.3 基于物理的方法部單元形狀度量低于一定的限值才轉(zhuǎn)到優(yōu)化的方法。這種改進(jìn)網(wǎng)格質(zhì)量的方法是通過(guò)模擬物理上任一節(jié)點(diǎn)和它周圍的節(jié)點(diǎn)之間存在著與距離有關(guān)的或排斥力的方法來(lái)重定位節(jié)點(diǎn)的。網(wǎng)格中或排斥力，重定位的結(jié)果是該點(diǎn)處于使它所受的力處于平衡的位置。Lohner 等96用一個(gè)相互作用的彈簧系統(tǒng)來(lái)模擬相鄰節(jié)點(diǎn)之間的關(guān)系，Shimada 等97和sen 等98把節(jié)點(diǎn)當(dāng)作氣泡的中心，重新定位這些氣泡的位置使它

46、們處于平衡，改變這些粒子（氣泡）之間力的大小和方向可以改變單元等向性和尺寸的大小。3.5.1.4 中間節(jié)點(diǎn)定位法多數(shù)平滑方法都是針對(duì)單元的角點(diǎn)進(jìn)行的。對(duì)存在中間節(jié)點(diǎn)的單元，Salem 等99提出了重定位中間節(jié)點(diǎn)來(lái)提高單元質(zhì)量的方法，這種方法計(jì)算環(huán)繞中間節(jié)點(diǎn)的一個(gè)稱為中間節(jié)點(diǎn)容許空間（Mid-node Admissible Space，MAS）的區(qū)域，在這個(gè)區(qū)域內(nèi)可以安全地移動(dòng)中間節(jié)點(diǎn)，保持或改進(jìn)單元的質(zhì)量。3.5.2 單元細(xì)分單元細(xì)分定義為在網(wǎng)格上完成的任何能有效減小單元尺寸的操作。減小單元尺寸可能是為了能得到更詳細(xì)的答數(shù)據(jù)，也可能僅僅是為了改善局部單元的質(zhì)量。有些細(xì)分方法本身可以認(rèn)為是網(wǎng)格生

47、成算法，從一個(gè)粗分的網(wǎng)格開始，用一個(gè)網(wǎng)格細(xì)分過(guò)程細(xì)分直到達(dá)到希望的節(jié)點(diǎn)密度。細(xì)分算法通常被作為自適應(yīng)分析的一部分，由先前的解答結(jié)果提供用于網(wǎng)格細(xì)分的準(zhǔn)則。圖3-15 局部四邊形網(wǎng)格細(xì)分Fig. 3-15 Exle of local quad refinement由于四邊形六面體網(wǎng)格本身的特性，通常不能應(yīng)用點(diǎn)和邊等分方法來(lái)細(xì)分網(wǎng)格。四邊33形和六面體網(wǎng)格細(xì)分的主要方法是基于一套預(yù)定義的模板來(lái)分解單元。文獻(xiàn)100描述了一系列的單元分解模板和算法，圖3-15顯示了一個(gè)四邊形局部細(xì)分的例子。為了維持網(wǎng)格的一致性，某些四邊形和六面體細(xì)分方法場(chǎng)需要引入三角形單元或四面體、五面體等單元。3.6 帶寬優(yōu)化有限

48、元分析等許多科學(xué)和工程的數(shù)值計(jì)算問(wèn)題都需要求解一個(gè)線性代數(shù)方Kx f(3-13)這里已知向量 f 和方陣 K，求解得到未知的向量 x。在有限元分析中，K 是一個(gè)稀疏矩陣（矩陣中的零元素遠(yuǎn)多于非零元素），矩陣 K 中非零元素的位置取決于未知量排列的次序，在有限元分析中未知量排列的次序取決于節(jié)點(diǎn)的次序。有限元分析中導(dǎo)出的剛度矩陣等的階數(shù)很高，一般是對(duì)稱矩陣，完全需要很大的空間，一般采用帶寬技術(shù)來(lái)剛度矩陣。帶寬技術(shù)的量取決于矩陣中主對(duì)角元素與該行非零量的最大距離（按行為 1 時(shí)，對(duì)節(jié)點(diǎn)進(jìn)行重新），而帶寬與單元節(jié)點(diǎn)的最大差值直接相關(guān)。當(dāng)每個(gè)節(jié)點(diǎn)的度等價(jià)于對(duì)方程(3-13)進(jìn)行重新排序；當(dāng)每個(gè)節(jié)點(diǎn)的度為

49、 m 時(shí)，則每個(gè)節(jié)點(diǎn)對(duì)應(yīng)于 m 個(gè)方程，對(duì)節(jié)點(diǎn)重新小帶寬的結(jié)果是一樣的。只是一個(gè)比對(duì)度重標(biāo)號(hào)小 O(m)，其減減小帶寬的算法主要有三種：直接減小算法（Direct Reduction Algorithms）、圖理論算法（Graph Theory Algorithms）和混合算法（Hybrid Algorithms）。直接減小算法是通過(guò)直接交換矩陣中的行和列來(lái)減小帶寬的，這種方法可以得到好的結(jié)果，但計(jì)算花費(fèi)的時(shí)間很大，這種算法的例子可參考文獻(xiàn)101。圖理論算法和圖理論算法與直接減小算法的混合方法計(jì)算效率高、結(jié)果好，以下主要對(duì)后兩種方法進(jìn)行說(shuō)明102,103,104,105,106,107,108

50、,109。3.6.1 有關(guān)的定義最合適的方程序號(hào)重標(biāo)算法取決于求解方程所使用的求解方法。有限元分析中根據(jù)剛度矩陣等所使用的方式，可分為二維等帶寬方法、一維變帶寬方法和。為評(píng)估算法的有效性，需要定義和矩陣有關(guān)的參數(shù)。設(shè)矩陣 K 為 n 階對(duì)稱方陣，定義第 i 行的行帶寬i(K)為該行對(duì)角線上的元素到該行最后一個(gè)不為零的元素的列數(shù)（包括對(duì)角元素） K maxj i 1, j i | A 0(3-14)iij這里假定 Aii 0，0 i n。矩陣 K 的帶寬(K)（Bandwidth）和廓 p(K)（Profile）定義為 K max(3-15)npK i K i1(3-16)34在矩陣 K 第 i

51、 行中，如果第 j 列滿足 j i，同時(shí)在任一行 k i 中存在 Kkj 0，則說(shuō) i 行中的第j 列是活動(dòng)的。用 ci(K)表示行 i 中的活動(dòng)列數(shù)，則矩陣 K 的最大波前定義為wmaxK maxci K , 0 i n方陣 K 是對(duì)稱矩陣，故有(3-17)nnpK i K ci K (3-18)i 1i 1平均波前 wavg(K)和均波前 wrms(K)定義為pK 1nwK c K i 1(3-19)avginnn1nwK c K i2(3-20)rmsi 1從以上的定義中，有wavg wrms wmax n(3-21)定義節(jié)點(diǎn) i 的度 di 為連接到節(jié)點(diǎn) i 的節(jié)點(diǎn)數(shù)目，如果不計(jì)節(jié)點(diǎn)的

52、度，di 是矩陣 K 中的 i行中除對(duì)角元素以外的非零項(xiàng)數(shù)。這里暗含著在四邊形或六面體單元中處于對(duì)角的節(jié)點(diǎn)是相互連接的。最大的節(jié)點(diǎn)度為m maxdi , 1 i n獨(dú)立的邊數(shù) e 定義為對(duì)角線以上非零元素的個(gè)數(shù)(3-22)1ne d(3-23)i2i1這樣矩陣 K 中總的非零元素的個(gè)數(shù)為 2en，矩陣 K 的密度 2e n(3-24)n 2需要注意的是在i 和 ci 的定義中包含矩陣 K 的對(duì)角元素。根據(jù)這些量的定義，可以方便的將各種參數(shù)從一種表示（包括對(duì)角元素）轉(zhuǎn)換為另一種表示（不包括對(duì)角元素）。同樣地，這里所說(shuō)的矩陣 K 的階 n 有時(shí)就是指節(jié)點(diǎn)的數(shù)目。在一般的有限元分析中，每個(gè)節(jié)點(diǎn)有數(shù)個(gè)

53、度，例如如果每個(gè)節(jié)點(diǎn)有 6 個(gè)度，則在沒(méi)有約束的情況下，wavg，wrms，wmax是僅考慮網(wǎng)格節(jié)點(diǎn)時(shí)對(duì)應(yīng)的數(shù)值的 6 倍。353.6.2 來(lái)源于圖理論的有關(guān)定義定義一個(gè)圖 G = (V, E)，V = v1, v2, , v n是頂點(diǎn)(Vertex)的有限非空集，vi, vjE，i j 定義了稱為邊(Edge)的有限集，V = V(G)，E = E(G)。這里 G 是一個(gè)無(wú)向圖（Undirected Graph）。如果vi, vj，則稱頂點(diǎn) vi 和 vj 是相鄰接的。對(duì)于有限元分析來(lái)說(shuō)，V 相當(dāng)于節(jié)點(diǎn)集合，E 相當(dāng)于節(jié)點(diǎn)之間的聯(lián)系。圖 G(V, E)的子圖（Subgraph）S(V, E)

54、定義為 SG，VV，EE。一個(gè)頂點(diǎn) v 的度（Degree）g(v)是與 v 相連的邊數(shù)。一個(gè)長(zhǎng)度為 t 的路徑（Path）是邊的序列v1,v2，v2,v3， vt-1,vt。兩個(gè)頂點(diǎn) u 和 v 之間存在一路徑，則說(shuō) u 和 v 是連通的（Connected）。u，v 之間的距離（Distance） (u, v)是連通 u，v 兩頂點(diǎn)的最短路徑。如果圖 G 中的任一頂點(diǎn)對(duì)是連通的，則圖 G 是連通的，否則是不連通的。如果圖 G 有 n 個(gè)頂點(diǎn) V = v1, v2, , vn，一對(duì)一的f 將v1, v2, , vn為1，2，n，f 稱為圖 G 的標(biāo)號(hào)（Numbering）。對(duì)每一個(gè)標(biāo)號(hào) f，

55、定義相對(duì)于 f 的圖 G 的帶寬為 f G max f v1 f v2 , v1, v2 E(3-25)一個(gè)圖 G 的層次結(jié)構(gòu)（Level Structure）L(G)是將圖 G 按如下方式分為子集 L1，L2， Ld與層 L1 中頂點(diǎn)相鄰接的每個(gè)頂點(diǎn)處于層 L1 或 L2 中；與層 Lj（1 j 0，做 f-i 步。選擇將被標(biāo)號(hào)的頂點(diǎn)。搜索優(yōu)先隊(duì)列，根據(jù) pt = maxpj，jQ定位頂點(diǎn) t，m 為頂點(diǎn) t的索引，qm = t。更新隊(duì)列和優(yōu)先級(jí)。通過(guò)設(shè)置 qm = qn 將頂點(diǎn) t 從優(yōu)先隊(duì)列中刪除，nn 1。如果頂點(diǎn)t 不是預(yù)活動(dòng)的，則進(jìn)行下一步；否則檢查與頂點(diǎn) t 鄰接的每個(gè)頂點(diǎn) x

56、并按照 px = px + W2增加它的優(yōu)先級(jí)（對(duì)應(yīng)于將頂點(diǎn) x 的度減小一個(gè)）。如果頂點(diǎn) x 是不活動(dòng)的，通過(guò)設(shè)置 nn + 1，qn = x 以預(yù)活動(dòng)狀態(tài)將 x優(yōu)先隊(duì)列中。h.標(biāo)號(hào)下一個(gè)頂點(diǎn)。置頂點(diǎn)計(jì)數(shù) ii + 1，用 i 標(biāo)號(hào)頂點(diǎn) t；置 rt = i，R 為新頂點(diǎn)標(biāo)號(hào)表。置頂點(diǎn) t 的狀態(tài)為過(guò)活動(dòng)的。i. 更新優(yōu)先級(jí)和隊(duì)列。檢查與頂點(diǎn) t 鄰接的每個(gè)頂點(diǎn) x，如果 x 不是預(yù)活動(dòng)的，置為不活動(dòng)；否則，置頂點(diǎn) x 為活動(dòng)狀態(tài)，px = px + W2，檢查與頂點(diǎn) x 相鄰接的每個(gè)頂點(diǎn) y。如果頂點(diǎn) y 不是過(guò)活動(dòng)的，增加它的優(yōu)先級(jí) py = py + W2；如果 y 是不活動(dòng)的，通過(guò)設(shè)

57、置 nn +1，qn = y 以預(yù)活動(dòng)狀態(tài)將 y優(yōu)先隊(duì)列中。j. 結(jié)束。新的頂點(diǎn)標(biāo)號(hào)保存在 R 中，rt 是頂點(diǎn) t 標(biāo)號(hào)。上面的算法中權(quán) W1 和 W2 反映了當(dāng)前頂點(diǎn)與端部頂點(diǎn)的距離和當(dāng)前頂點(diǎn)度的相對(duì)重要性，根據(jù)經(jīng)驗(yàn)一般取 W1= 1 和 W2 = 2，并且 W2 W1。以上的重標(biāo)號(hào)算法適用于使用帶寬求解的求解方案。使用求解時(shí)，由于方程的集合和刪除是以單元依次進(jìn)行的，故必須決定一個(gè)高效的單元次序。按已優(yōu)化的節(jié)點(diǎn)，以單元中最小節(jié)點(diǎn)的升序分類和標(biāo)號(hào)單元，可以導(dǎo)出一個(gè)好的單元標(biāo)號(hào)。這樣做保證了刪除方程是以類似于優(yōu)化的節(jié)點(diǎn)的次序進(jìn)行的，實(shí)踐證明也是有效的。Sloan 算法的效果較好，且適應(yīng)帶寬解法

58、和波前解法，編程也較容易。3.6.3.4 混合方法Armstrong104混合算法先使用 GPS 算法將節(jié)點(diǎn)重標(biāo)號(hào)，然后使用行列交換的的方法來(lái)改進(jìn)標(biāo)號(hào)的次序。Armstrong 觀察到基于圖理論的算法能快速給出一個(gè)粗略的優(yōu)化節(jié)點(diǎn)標(biāo)號(hào)，但不能象行列交換算法那樣得到更好的、詳細(xì)的節(jié)點(diǎn)。Armstrong 的混合算法可以得到好于 GPS算法的結(jié)果。3.7 小結(jié)本章對(duì)有限元網(wǎng)格生成技術(shù)，特別是三維非結(jié)構(gòu)六面體網(wǎng)格生成技術(shù)進(jìn)行了全面的總結(jié)，并對(duì)有關(guān)的技術(shù)進(jìn)行了較為詳細(xì)的描述。對(duì)有限元網(wǎng)格后處理所使用的技術(shù)，如網(wǎng)格節(jié)點(diǎn)位置平滑、帶寬優(yōu)化等技術(shù)進(jìn)行了總結(jié)，并詳細(xì)的描述了有關(guān)的算法。40K. Ho-Le, F

59、inite Element Mesh Generation Methods: A Review and Classification, Computer Aided Design, 1988, 20:2738.M. S. Shepherd, K. R. Grice, J. O. Lo and W. J. Schroeder, Trends in Automatic ThreeDimenal Mesh Generation, Computers & Structures, 30(1-2):42142924 D. A. Field, The Legacy of Automatic Mesh Gen

60、eration from Solid Ming, Computer Aided Geometric Design, 1995, 12:65167325，有限元網(wǎng)格生成方法發(fā)展綜述，計(jì)算機(jī)輔助設(shè)計(jì)和圖形學(xué)學(xué)報(bào)，1997，9(4):37838326，(2):2831，有限元網(wǎng)格自動(dòng)生成典型方法及發(fā)展方向，計(jì)算機(jī)輔助設(shè)計(jì)和制造，1996，27 閔，有限元網(wǎng)格劃分技術(shù)，計(jì)算機(jī)研究與發(fā)展，1995，32(7):3742，28 C. Gay, P. Montmitonnet, N. Soyris and J. L. Chenot, Validation of a 3-D Finite Ele ment P