1、第14卷第4期2007年8月塑性工程學報J OU RNAL OF PL ASTICIT Y EN GIN EERIN GVol 114No 14Aug 12007全自動四邊形有限元網格生成-I :區(qū)域分解法3(山東大學模具工程技術研究中心,濟南250061馬新武1趙國群3(山東建筑大學機電工程學院,濟南250101王芳2摘要:提出了一種適用于金屬體積成形有限元分析的全自動四邊形網格生成算法區(qū)域分解法,對該算法的一些關鍵技術,如幾何形狀的表示、關鍵點的確定、網格密度的生成、節(jié)點的生成、邊界單元的生成、剖分準則、模式匹配、光滑處理以及半帶寬的優(yōu)化等進行了詳細的介紹。以該算法為基礎,成功開發(fā)了全自動

2、四邊形網格生成軟件AU TOM ESH 。并給出了一些網格劃分的實例,驗證該算法的可靠性。關鍵詞:區(qū)域分解法;網格生成;四邊形網格;密度控制中圖分類號:TB115文獻標識碼:A 文章編號:100722012(200704201052053山東省自然科學基金資助項目(Q2004F01;國家自然科學基金資助項目(50375087。馬新武E 2mail :maxinwu sdu 1edu 1cn作者簡介:馬新武,男,1975年生,副教授,主要從事材料成形過程數值模擬、優(yōu)化設計及軟件開發(fā)研究收稿日期:2006203216;修訂日期:2007204221前言盡管塑性有限元理論和技術已有很大的發(fā)展,但由于

3、塑性成形自身的特點,使得有限元法在實際應用中仍存在許多具體問題和數學處理上的困難。塑性成形的特點是變形量大,模具表面和變形材料之間的相對運動也非常大。因此,利用有限元法模擬塑性成形過程時,有限元網格很容易發(fā)生單元畸變或與模具干涉,導致有限元分析過程中斷。必須定義新的網格系統(tǒng),才能重新開始有限元分析過程。實際上,沒有網格重劃分,也就不可能模擬塑性成形過程。有限元網格自動劃分及重劃分技術已是影響有限元模擬系統(tǒng)通用化的“瓶頸”之一。目前影響國內體積成形有限元軟件推廣的主要原因之一,就是缺少可靠、高效、適合體積成形特點的網格自動生成算法和程序。生成四邊形網格的方法很多,其中最有代表性和使用最多的當屬幾

4、何分解法。在幾何分解法中Loop 15和Paving 6方法是生成單元質量較好的兩種方法,但這兩種方法的缺點是很難控制網格密度,而體積成形的特點是要求網格密度易于控制:在工件與模具接觸的部位要求網格密度大些,單元細小,以使工件邊界和模具貼合;而在幾何形狀較平直的部位網格密度可以小些,且網格疏密應過渡均勻。本文在開發(fā)體積成形有限元軟件CASFORM 7的過程中,針對體積成形的特點,提出了一種密度控制和幾何分解相結合的四邊形網格自動生成方法,稱為區(qū)域分解法;并以此為基礎,成功地開發(fā)出了四邊形網格自動生成程序AU TOM ESH 8。本文中主要介紹區(qū)域分解法及其關鍵技術,下一篇中介紹網格密度控制和生

5、成方法。1區(qū)域分解法在區(qū)域分解法中,首先根據物體的幾何邊界曲率大小、網格密度過渡大小要求以及要劃分單元的數目,在整個物體上生成過渡均勻的網格。然后根據密度分布,在物體上生成四邊形網格。網格的生成過程可以用圖1所示的簡單例子加以說明。首先圖1網格生成過程示意圖Fig 11Schematical illustration of the meshgeneration process根據密度分布,在物體的邊界上生成邊界節(jié)點,通過連接兩個邊界節(jié)點,把這個區(qū)域剖分成兩個子區(qū)域。再根據密度分布在剖分線上增加新的節(jié)點。以遞歸的方式對每個子區(qū)域進行剖分,直到所有的子 區(qū)域不可再分為止,即每個子區(qū)域包含6個或4個

6、節(jié)點。在剖分完成后,每一個6節(jié)點子區(qū)域映射到3種基本模式之一上,生成最后的單元。需要注意的是,為了生成四邊形網格,邊界節(jié)點總數必須是偶數。因此,不僅在生成初始的邊界節(jié)點時,而且在剖分線上生成節(jié)點時都要滿足這個要求。這樣區(qū)域和子區(qū)域都可以用邊界節(jié)點環(huán)來表示。區(qū)域分解法用到的主要關鍵技術包括:幾何形狀的表示、關鍵點的確定、網格密度的生成、節(jié)點的生成、邊界單元的生成、剖分準則、模式匹配、光滑處理以及半帶寬的優(yōu)化等。2關鍵技術211幾何形狀的表示考慮到大多數模具和工件形狀都可以由直線和圓弧來表示,因此,對于物體的幾何形狀,本文采用節(jié)點坐標和節(jié)點處的過渡圓角半徑來表示,并且規(guī)定兩條直線段交點的過渡圓角半

7、徑為0,這樣既可以用于初始的幾何形狀輸入,也可用于再劃分時的工件邊界形狀的表示。對于圓弧部分還需要根據精度要求用直線段來擬合,因此幾何形狀最終是由首尾相連的直線段組成的多邊形來表示。對于復雜形狀曲線,例如B 樣條曲線,可以在精度許可的前提下,用直線段來逼近。212關鍵點的確定表示物體形狀的多邊形中能體現物體形狀特征的關鍵節(jié)點,例如尖角處節(jié)點,稱為關鍵點。確定關鍵點后,物體形狀就可以用含有最少節(jié)點的特征多邊形來表示。本文采用一種叫做多邊形逼近的算法9來確定特征多邊形的關鍵點。213網格密度的生成有限元自適應技術要求網格劃分程序能夠根據物體的幾何形狀和場變量的分布,合理布置網格密度。從直觀上看,網

8、格密度大的地方,網格單元尺寸較小;網格密度小的地方,網格單元尺寸較大。從數學上可以定義網格密度為網格單元長度的倒數。本文采用Laplace 方程來表示物體上的網格密度分布函數。網格密度的生成過程是:首先在物體上劃分均勻的柵格,根據網格密度的要求,在物體邊界上和內部指定網格密度,采用有限元求解技術,在整個物體上生成過渡均勻的網格密度。這樣物體上任一點的密度都可通過柵格節(jié)點上的密度插值得到。關于網格密度生成和控制的詳細討論,請見本文的下篇網格密度控制。214節(jié)點的生成生成網格密度之后,就可以根據網格密度分布在物體邊界上生成節(jié)點。節(jié)點的生成過程是在物體特征多邊形每條邊上,即兩個關鍵點之間生成節(jié)點,然

9、后按生成的順序依次給節(jié)點編號。在兩個關鍵點之間生成節(jié)點需要兩個步驟。第一步需要計算生成節(jié)點的數目,第二步根據網格密度分布決定各個節(jié)點的位置。要計算生成節(jié)點的數目,需要沿著邊界在兩個關鍵點之間對網格密度積分,即Lu (s d s =R N -1(1式中R 積分結果N 包括兩個端節(jié)點在內生成的節(jié)點數L 兩個關鍵點間的連線u (s 邊界上的網格密度d s 沿邊界的微分在實際數值計算時,式(1可以按下面的方法進行計算:將積分連線L 分成若干等長的小線段,將每一小線段的中點處密度值作為該段的平均密度值,用每一小段的平均密度值乘以該段的長度,然后對這些乘積求和,就可以得到近似的積分結果R 。積分結果R 為

10、一實數,需要四舍五入到最接近的整數N -1。另外,上述已提及到邊界節(jié)點總數必須是偶數個,所以為了滿足這個要求,必要時還需要對N 值進行調整。正是由于存在舍入誤差,所以最后生成的網格單元數目和設定值之間存在一定的偏差。確定了要生成節(jié)點的數目N 后,下一步就要確定這些節(jié)點的位置。除了第1個和第N 個節(jié)點的位置固定之外,其余的節(jié)點位置都是未知的?？梢詫⒐?jié)點的位置(相對起點的長度作為一維變量,則該變量應滿足控制方程u (s d s =A d=C (2式中u (s 沿邊界的密度函數d s 實際空間上的微分d 計算空間上的微分(在計算空間上,節(jié)點是均勻分布的A 比例常數C 常數式(2可寫成兩階微分形式為d

11、 du (s d sd =0(3601塑性工程學報第14卷對于式(3所示的微分方程,可以采用有限元迭代技術求得節(jié)點位置的數值解。將計算空間平均分成N -1個單元,共N 個節(jié)點。經過單元分析,可得到單元剛度方程為uS i +S i+ 121- 1-11S i S i+1=00(4式中S i 第i 個節(jié)點的位置S i +1第i +1個節(jié)點的位置對所有單元的剛度方程進行組裝可得到整體剛度方程,將第1個節(jié)點位置S 1=0,第N 個節(jié)點位置S N =L (兩個關鍵點的長度作為邊界條件,其余節(jié)點的初始位置假定是均勻分布的,對整體剛度方程進行迭代求解,當求解過程收斂時就可得到各節(jié)點的位置。圖2a 為一網格密

12、度分布函數 ,利用上述方法計算得到生成的節(jié)點數目N =37,各節(jié)點的位置如圖2b 所示。圖2根據網格密度布置及節(jié)點位置a 網格密度分布;b 節(jié)點的位置Fig 12Nodes located based on mesh density215邊界單元的生成由于塑性成形的特點,用有限元法模擬工件成形時,邊界單元最容易發(fā)生變形,網格發(fā)生畸變,導致網格重劃分。如果在網格劃分時工件邊界單元網格質量較高,則明顯地提高整個網格系統(tǒng)的抗畸變能力,減少模擬過程中的網格重劃分次數,提高模擬精度。在網格生成方法中,Paving 方法將要劃分的區(qū)域由邊界逐層向內部推進生成網格單元,邊界網格單元形狀規(guī)整,質量較高。本文借

13、鑒Paving 方法的思路,在生成所有的邊界節(jié)點之后,通過節(jié)點偏移的方法,在邊界生成一層網格單元。這樣由偏移的節(jié)點組成的節(jié)點環(huán)在物體的內部又形成一個新的區(qū)域。下一步的任務就是要將這個區(qū)域剖分成兩個子區(qū)域。216剖分準則從區(qū)域邊界節(jié)點環(huán)中選擇兩個節(jié)點,通過連接這兩個節(jié)點形成一條剖分線就可以將區(qū)域一分為二。如何選擇剖分線不僅影響到最后生成的網格質量,也影響到網格生成的效率。目前大部分文獻都是以剖分角度、剖分線的長度,以及剖分后兩子區(qū)域的面積比作為選擇剖分線的準則。如圖3所示,連接節(jié)點i ,j 生成一條剖分線i j ,將區(qū)域分成兩個子區(qū)域,同時在節(jié)點i ,j 處形成4個剖分角ij 1,ij 2和ji

14、 1,ji 2。原則上4個剖分角越接近直角,剖分線的長度越短,兩子區(qū)域面積比越接近1,那么這樣的剖分線就是最佳的剖分線。通常的做法是將上述3個標準進行量化,采用加權平均的方法得到一個量化值,對所有可能的剖分線比較其量化值,選擇量化值最小的剖分線作為最佳的剖分線。實踐中,上述方法存在如下的缺點:13個標準量化權值的最優(yōu)值很難確定,文獻中給出的權值各不相同,所給出的經驗值,并不一定適用所有的情況;2計算時間長,對于有n 個節(jié)點的節(jié)點環(huán),可能的剖分線條數為n (n -3/2,由于對每條可能的剖分線都要計算剖分角、剖分線長度以及子區(qū)域的面積,使得計算時間特別長,影響了網格劃分的效率。本文提出一種確定最

15、佳剖分線的新方法,在滿足網格劃分質量的前提下,顯著地提高網格劃分的效率。該方法確定最佳剖分線的步驟是:1如果節(jié)點i 處的角度i <60°,該節(jié)點就不能與其他的節(jié)點連接生成剖分線;2如果連接節(jié)點i ,j 形成的4個剖分角ij 1,ij 2,ji 1,ji 2中有一個小于30°,節(jié)點i ,j 之間就不能生成剖分線;3如果連接節(jié)點i ,j 的線段與邊界相交,則節(jié)點i ,j 之間就不能生成剖分線;4除了前面幾種情況外,連接節(jié)點i ,j 就可以確定一條剖分線。節(jié)點i ,j 配對后就不再與其他節(jié)點配對確定剖分線。計算剖分線i j 的長度和剖分后兩子區(qū)域的面積比(1;5繼續(xù)尋找其他

16、可能的剖分線,并計算各剖分線的長度和剖分后兩子區(qū)域的面積比;6對所有可能的剖分線按長度和子區(qū)域的面積比從小到大排序,選擇子區(qū)域面積比大于013,且長度最短的剖分線作為最佳的剖分線;7按照上述步驟無法確定剖分線時,將邊界節(jié)點環(huán)壓入棧中,用后面介紹的模式匹配方法進行處理。701第4期馬新武等:全自動四邊形有限元網格生成-I :區(qū)域分解法 確定最佳剖分線后,就可以按照網格密度要求在剖分線上生成新節(jié)點,于是原區(qū)域一分為二,形成兩個新的節(jié)點環(huán)。對于每個子區(qū)域,按照上述方法進行處理,直到所有子區(qū)域不可再分為止 。圖3剖分線Fig 13Splitting line217模式匹配當物體的幾何形狀分解完成之后,

17、生成的六節(jié)點區(qū)域可以映射到3種基本模式之一,生成最后的四邊形單元,對于那些多于六節(jié)點而又不滿足剖分準則的區(qū)域,可以通過遞歸的方式將該區(qū)域分解成為3種基本模式。圖4為六節(jié)點區(qū)域的3種基本模式,對于模式1連接相對邊上的兩個節(jié)點,將區(qū)域分解成兩個四邊形單元;對于模式2和3,分別在區(qū)域內部增加1個和2個新節(jié)點,將區(qū)域分解成3個和4個四邊形單元。圖5為多于六節(jié)點的區(qū)域分解成基本模式的示意圖 。圖43種基本模式Fig 14The three basic patterns圖5模式匹配過程示意圖Fig 15Schematic illustration of patternmatching procedure2

18、18光滑處理和半帶寬優(yōu)化由于直線剖分和模式匹配等原因,區(qū)域分解法生成的最后單元形狀可能較差。為了改進單元形狀,本文用到了Laplace 光滑算法對最后生成的網格進行光滑處理。同樣,用該方法生成的網格進行有限元分析時,其剛度矩陣的帶寬分布較差。為獲得計算效率高的網格,必須減少網格的帶寬,對節(jié)點編號進行優(yōu)化。3網格劃分實例基于區(qū)域分解法開發(fā)了全自動四邊形網格生成軟件AU TOM ESH 。應用AU TOM ESH 軟件,進行了大量的實例劃分,并結合體積成形仿真軟件CASFORM 進行了大量的模擬分析,表明了該方法的可靠性。圖6為模擬一復雜擠壓過程時采用區(qū)域分解法對工件進行網格再劃分的實例 。圖6網

19、格劃分實例Fig 16Examples of mesh generation4網格劃分質量和效率反映平面四邊形網格質量的指標通常有兩個:801塑性工程學報第14卷一是網格單元相鄰邊的長度比,一是網格單元各內角的角度。對于均勻劃分的網格,單元鄰邊長度比越接近1,單元各內角越接近90°,網格質量就越高;對于密度變化的網格,要求網格尺寸大小從密度大的區(qū)域到密度小的區(qū)域應過渡均勻。考慮到要劃分網格的區(qū)域形狀復雜程度不一樣,通常要求網格單元的內角在30°150°之間。從大量的劃分實例結果來看,采用區(qū)域分解法生成的網格單元的內角通常都在30°150°之間,

20、尤其是邊界單元,其質量更高。由于采用了密度控制技術,該方法生成的網格大小疏密過渡均勻。區(qū)域分解法生成網格的效率非常高,通常只需要輸入幾何形狀和要求劃分的網格數目及最大單元和最小單元的面積比兩個參數,即可自動生成過渡均勻的網格密度;并且根據網格密度自動生成網格,整個過程不需要人工的干預,這對于通常需要多次網格重劃的體積成形有限元模擬來說特別重要。區(qū)域分解法生成網格的速度也非?？?在Pentium III 微機上生成1000個單元僅需1s 3s 的時間。5結論本文提出了一種用于金屬體積成形有限元分析的四邊形網格生成算法區(qū)域分解法,該算法具有算法可靠、劃分速度快、劃分質量高、網格密度和單元數目易于控

21、制、輸入參數少、易于實現等優(yōu)點。區(qū)域分解法不僅可用來生成用于有限元分析的初始網格,也可用于有限元分析過程中網格畸變后的再劃分,因此特別適合于金屬體積成形有限元網格劃分。同時適用于其他結構有限元分析的網格生成。參考文獻1J A Talbert ,A R Parkinson.Development of an Auto 2matic ,Two 2dimension Finite Element Mesh Generator Using Quadrilateral Elements and Bezier Curve Bound 2ary Definition J .Int.J.Numer.Meth.

22、Engng ,1990.(29:155115672W T Wu ,S I Oh et al.Automated Mesh G eneration forForming Simulation I C .ASM E International Com 2puters in Engineering Conference Proceedings ,Boston ,Massachusetts ,1990.1,5073郭曉霞.特別適用于金屬成形模擬的四邊形網格劃分方法J .塑性工程學報,2003.(6:50544郭曉霞.四邊形有限元網格劃分方法二分法的改進J .塑性工程學報,2005.12(5:96100

23、5J Sarrate ,A Huerta.Efficient Unstructured Quadrilat 2eral Mesh GenerationJ .Int.J.Numer.Meth.Engng ,2000.(49:132713506T D Blacker ,M B Stephenson.Paving :A New Ap 2proach to Automated Quadrilateral Mesh Generation J .Int.J.Numer.Meth.Engng.,1991.(32:8118477馬新武,趙國群.體積成形有限元模擬軟件CASFORM的開發(fā)研究J .鍛壓技術,20

24、03.(1:9138趙國群,馬新武.全自動自適應四邊形有限元網格生成/再生成方法與軟件AU TOM ESH 22D/PCJ .金屬成形工藝,2003.(2:42449T Pavlidis.Algorithms for Graphics and Image Pro 2cessingM .Computer Science Press ,1982The autom atic generation of qu adrilateral mesh 2I :the method of dom ain divisionMA Xin 2wu 1ZHAO Guo 2qun 3(Die &Mold Engineering Research Center ,Shandong University ,Jinan ,Shandong 250061China WAN G Fang 2(College of Mechanical and Electronic Engineering ,Shandong Jianzhu University ,Jinan ,Sha