蒲軍平,姚宇龍
(浙江工業(yè)大學(xué) 建筑工程學(xué)院,浙江 杭州 310032)
平面有限元網(wǎng)格實(shí)用劃分方法研究
蒲軍平,姚宇龍
(浙江工業(yè)大學(xué) 建筑工程學(xué)院,浙江 杭州 310032)
網(wǎng)格的生成是平面有限元計(jì)算的前提和基礎(chǔ),但當(dāng)邊界不規(guī)則時(shí)容易生成劣質(zhì)單元,對(duì)計(jì)算造成影響.對(duì)比目前較為成熟的四邊形網(wǎng)格生成方法,開(kāi)發(fā)出以Visual Basic語(yǔ)言為基礎(chǔ)的結(jié)構(gòu)化有限元網(wǎng)格劃分程序,生成較為規(guī)則的四邊形單元.它具有操作簡(jiǎn)便,網(wǎng)格密度、單元數(shù)量可控制,生成的數(shù)據(jù)簡(jiǎn)單實(shí)用等特點(diǎn),特別適用于小型結(jié)構(gòu)的有限元分析.
網(wǎng)格生成;四邊形單元;網(wǎng)格密度;Visual Basic語(yǔ)言
有限元法是一種有效的工程分析方法,特點(diǎn)是可以對(duì)不規(guī)則邊界、材料各項(xiàng)異性、有裂縫或突變等復(fù)雜工況進(jìn)行較為精確的分析.由于它在解決實(shí)際問(wèn)題中簡(jiǎn)便有效,在工程和科技領(lǐng)域中已經(jīng)得到了廣泛的應(yīng)用.有限元網(wǎng)格生成是指如何在求解區(qū)域內(nèi)生成拓?fù)溆行Ш蛶缀斡行У臐M(mǎn)足有限元計(jì)算要求的網(wǎng)格,目前有限元網(wǎng)格的研究已取得了許多重要的成果,科研人員研究并總結(jié)出了一些成熟的算法在此基礎(chǔ)上研制了一批面向市場(chǎng)的商業(yè)工程化軟件.國(guó)內(nèi)的有限元網(wǎng)格技術(shù)起步晚于國(guó)外,但發(fā)展迅速,隨著計(jì)算機(jī)技術(shù)的發(fā)展,國(guó)內(nèi)的一些專(zhuān)家學(xué)者也提出了很多優(yōu)秀的網(wǎng)格生成的算法來(lái)解決工程實(shí)際中對(duì)有限元網(wǎng)格計(jì)算的要求[1-2].如何在現(xiàn)有的基礎(chǔ)上更好地實(shí)現(xiàn)有限元網(wǎng)格的前處理,除了對(duì)算法本身進(jìn)行研究,尋找精度更高,收斂更快,試用范圍更廣的理論外還要有與之相適應(yīng)的計(jì)算機(jī)程序匹配才能顯示出它的優(yōu)越性,二維有限元的前處理的核心就是進(jìn)行合理高效的網(wǎng)格剖分系統(tǒng)的設(shè)計(jì).筆者研究的目的在于找到一種適合處理小型復(fù)雜區(qū)域的平面網(wǎng)格劃分方法,進(jìn)而為之后的有限元計(jì)算做好前期準(zhǔn)備工作.
現(xiàn)有的有限元網(wǎng)格生成方法有數(shù)百種之多[3-5],但大致能歸納為兩大類(lèi),即結(jié)構(gòu)化有限元和非結(jié)構(gòu)化有限元網(wǎng)格.
結(jié)構(gòu)化有限元網(wǎng)格指網(wǎng)格內(nèi)部所有點(diǎn)的度均相同,具體說(shuō)就是內(nèi)部所有節(jié)點(diǎn)具有數(shù)量相等的毗鄰單元,結(jié)構(gòu)化有限元網(wǎng)格有許多優(yōu)點(diǎn):
1)邊界附近的單元與邊界可以對(duì)齊,對(duì)邊界的處理相對(duì)來(lái)說(shuō)比較簡(jiǎn)單,特別適用于應(yīng)力集中或者流體、磁場(chǎng)等對(duì)邊界條件敏感的有限元計(jì)算.
2)網(wǎng)格生成的數(shù)量、密度可人為控制且生成的單元質(zhì)量好.
3)數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)單、實(shí)用,特別適用于小型結(jié)構(gòu)的有限元計(jì)算.
當(dāng)然結(jié)構(gòu)化網(wǎng)格的缺點(diǎn)也是顯而易見(jiàn)的,一般來(lái)說(shuō)它只能對(duì)特定的單元進(jìn)行劃分,如果是復(fù)雜的區(qū)域,只能事先人為地將復(fù)雜區(qū)域劃分為若干個(gè)四邊形區(qū)域,然后對(duì)每個(gè)四邊形區(qū)域單獨(dú)劃分.通常情況下,四邊形單元是比較理想的二維結(jié)構(gòu)化網(wǎng)格,筆者所生成的即是結(jié)構(gòu)化有限元網(wǎng)格.
非結(jié)構(gòu)化有限元網(wǎng)格可以彌補(bǔ)結(jié)構(gòu)化網(wǎng)格難以解決的處理任意形狀和聯(lián)通區(qū)域的問(wèn)題,由于非結(jié)構(gòu)化網(wǎng)格生成技術(shù)比較復(fù)雜,算法和編程技術(shù)上的實(shí)現(xiàn)較為困難,目前采用比較成熟的非結(jié)構(gòu)化網(wǎng)格自動(dòng)生成的大多是三角形單元.
目前,平面有限元分析單元類(lèi)型有多種多樣,比較常見(jiàn)的是三角單元和四邊形單元[6-7],發(fā)展到今天三角形網(wǎng)格單元已基本成熟,基本上采用Delaunary三角算法就可得到滿(mǎn)意的結(jié)果.而四邊形網(wǎng)格的發(fā)展中雖然沒(méi)有一種比較統(tǒng)一的生成算法,但沿用四邊形單元進(jìn)行的有限元分析計(jì)算結(jié)果精度更高,所以現(xiàn)在的二維有限元分析方法都傾向于采用四邊形單元.算法的種類(lèi)多種多樣,根據(jù)其幾何特點(diǎn)可大致歸納為以下幾種:
前沿推進(jìn)法(Advancing front method)是較為常見(jiàn)的一種三角形網(wǎng)格生成方法,適用于結(jié)構(gòu)和非結(jié)構(gòu)化的三角形網(wǎng)格單元生成,經(jīng)過(guò)改良后也可用于四邊形網(wǎng)格單元.該方法的基本方法是,從區(qū)域的邊界開(kāi)始,按一定的方向形成一個(gè)前沿單元,并在此基礎(chǔ)上生成下一個(gè)滿(mǎn)足一定要求的單元,同時(shí)更新上一個(gè)前沿單元,以此重復(fù),直至充滿(mǎn)整個(gè)區(qū)域.使用前沿推進(jìn)法時(shí),可通過(guò)引入和改變單元尺寸函數(shù)來(lái)控制生成的子單元大小,在應(yīng)力集中或區(qū)域圖形有突變處加密單元,進(jìn)而讓網(wǎng)格劃分后的單元和數(shù)據(jù)更符合有限元計(jì)算的精度要求.
映射法是以手工的方式將復(fù)雜區(qū)域人為地劃為簡(jiǎn)單的區(qū)域,并將其映射到可以生成網(wǎng)格的規(guī)則區(qū)域.它的特點(diǎn)是算法簡(jiǎn)單、形成單元質(zhì)量高、速度快,既可生成結(jié)構(gòu)化網(wǎng)格亦可生成非結(jié)構(gòu)化網(wǎng)格.映射法將復(fù)雜的求解區(qū)域通過(guò)坐標(biāo)的轉(zhuǎn)化變成形狀規(guī)則的幾何圖形,然后在新的坐標(biāo)下對(duì)規(guī)則的區(qū)域進(jìn)行網(wǎng)格劃分,最后把網(wǎng)格劃分的結(jié)果轉(zhuǎn)化到原坐標(biāo)中,從而得到了經(jīng)過(guò)坐標(biāo)映射變化而來(lái)的網(wǎng)格劃分圖形.但是在曲面網(wǎng)格中產(chǎn)生的良好單元映射到原坐標(biāo)時(shí)會(huì)產(chǎn)生畸變.
柵格法一般可分為柵格疊合與四叉樹(shù)法,或兩者合在一起.先將柵格置放在物體上,節(jié)點(diǎn)既可布置在柵格規(guī)則點(diǎn)處,亦可在柵格單元內(nèi)隨機(jī)設(shè)置,可根據(jù)所需求解的區(qū)域調(diào)整,得到較為均勻的網(wǎng)格.再根據(jù)柵格與目標(biāo)圖形的求交以刪除外部的柵格,但保留內(nèi)部的柵格作為求解域的網(wǎng)格.刪除邊界外部柵格后,需調(diào)整節(jié)點(diǎn)的位置來(lái)適應(yīng)邊界條件的要求.
筆者所采用的Visual Basic程序語(yǔ)言,相對(duì)于C++,F(xiàn)ortran等計(jì)算機(jī)語(yǔ)言具有結(jié)構(gòu)簡(jiǎn)單,可視化強(qiáng)等優(yōu)勢(shì).由于傳統(tǒng)的可視化界面大都采用輸入坐標(biāo)等方式形成求解區(qū)域,在形成復(fù)雜圖形時(shí)非常不方便,例如一個(gè)簡(jiǎn)單的四邊形便要輸入8個(gè)坐標(biāo).為了彌補(bǔ)這一不足,筆者采用了VB交互式CAD技術(shù)[8].交互式CAD技術(shù)比一般的計(jì)算機(jī)繪圖技術(shù)要復(fù)雜一些,在交互式CAD系統(tǒng)中,不僅可以實(shí)現(xiàn)用鼠標(biāo)繪圖,還可以選擇單元,對(duì)選中的單元進(jìn)行復(fù)制、移動(dòng)、刪除、劃分等操作,為了實(shí)現(xiàn)這樣一些功能,交互式CAD系統(tǒng)不僅要考慮圖形的數(shù)學(xué)模型,還要在數(shù)據(jù)結(jié)構(gòu)、內(nèi)存管理等方面做出整體籌劃.傳統(tǒng)的數(shù)據(jù)輸入方式是以文件的形式將圖形控制點(diǎn)的坐標(biāo)輸入,再按一定的順序連接起控制點(diǎn).在這個(gè)過(guò)程中,具體的繪圖過(guò)程無(wú)法控制,是黑箱,只有輸入和輸出是可見(jiàn)的,為了得到想要的求解區(qū)域必須在稿紙上對(duì)每個(gè)控制點(diǎn)的坐標(biāo)進(jìn)行精確的求解.當(dāng)求解的區(qū)域較為復(fù)雜,圖形的生成就顯得異常繁瑣.另外,這種繪圖方式得到的圖形不能進(jìn)行修改和編輯,必須要修改數(shù)據(jù)文件.而這些需要修改的數(shù)據(jù)又很零碎的存放在大的數(shù)據(jù)文件中,很難找出相對(duì)應(yīng)的圖元數(shù)據(jù).在傳統(tǒng)的有限元計(jì)算中經(jīng)常會(huì)遇到數(shù)據(jù)的輸入和修改帶來(lái)的一系列問(wèn)題,特別是三維的有限元計(jì)算,要在三維空間中找到出現(xiàn)問(wèn)題的點(diǎn),然后查找相應(yīng)的數(shù)據(jù)文件,對(duì)其修改.交互式CAD的繪圖和建模技術(shù)對(duì)數(shù)據(jù)輸入和修改帶來(lái)了突破性的進(jìn)展,筆者采用的交互式有限元網(wǎng)格生成,可以在屏幕上移動(dòng)鼠標(biāo)光標(biāo),根據(jù)光標(biāo)處的坐標(biāo)進(jìn)行繪圖,這樣整個(gè)繪圖的過(guò)程就可以方便地控制.由于采用了特殊的數(shù)據(jù)結(jié)構(gòu)保存數(shù)據(jù),在圖形的形成和劃分過(guò)程中可以方便地對(duì)圖元進(jìn)行選擇、修改和編輯.
在網(wǎng)格劃分時(shí)采用VB中的“類(lèi)”這一特殊的形式作為描述對(duì)象,可以有效地集成各種具有類(lèi)似屬性和方法的代碼,并將其存在內(nèi)存中,方便下次調(diào)用時(shí)找出對(duì)應(yīng)的代碼并加以改動(dòng).例如,用鼠標(biāo)繪制直線(xiàn)和圓弧時(shí)同樣是單擊左鍵,但出現(xiàn)的效果分別是一條直線(xiàn)和一段圓弧,兩者可以放在同一個(gè)“CComand”類(lèi)里面,而在直線(xiàn)“CCreateLine”和圓弧“CCreateCircle”兩個(gè)類(lèi)里分別添加不同的代碼.使用這樣的方式可有效地規(guī)避VB不提供對(duì)象繼承的缺點(diǎn),讓一些類(lèi)似的模塊公用,有效地節(jié)省了代碼空間.
我們的生成思路是先將一個(gè)復(fù)雜的求解區(qū)域人為地劃分為若干個(gè)四邊形或三角形區(qū)域,再對(duì)每個(gè)子單元進(jìn)行網(wǎng)格劃分.由于四邊形網(wǎng)格在有限元計(jì)算精度上的優(yōu)勢(shì),所以采用了四邊形單元作為基本單元,而將少數(shù)三角形單元視為一邊退化的特殊四邊形.要提高有限元模擬的精度就必須提高網(wǎng)格劃分的質(zhì)量,譬如在構(gòu)件薄弱處、易產(chǎn)生應(yīng)力集中的部位對(duì)網(wǎng)格進(jìn)行加密,即求得合理的網(wǎng)格密度分布.由于要考慮到幾何形狀、劃分的單元數(shù)以及網(wǎng)格密度的過(guò)渡等諸多情況,一般可采用Laplace方程作為劃分區(qū)域內(nèi)網(wǎng)格密度函數(shù)的指導(dǎo),使得網(wǎng)格密度分布連續(xù)且過(guò)渡光滑.
其中:u為網(wǎng)格密度分布函數(shù);Ω為計(jì)算區(qū)域;Γ為計(jì)算區(qū)域的邊界;n為邊界Γ的法向方向.
由于生成的四邊形單元是獨(dú)立的單元,而每個(gè)獨(dú)立單元的網(wǎng)格劃分算法又略有差異,導(dǎo)致它們之間邊界的處理要格外謹(jǐn)慎.計(jì)算機(jī)可能會(huì)把邊界上非常接近的點(diǎn)視為兩個(gè)點(diǎn)從而令單元和節(jié)點(diǎn)號(hào)的生成產(chǎn)生混亂,處理的方法可將相近的兩個(gè)網(wǎng)格點(diǎn)做x,y向比較,小于某一界限值的即認(rèn)為是一個(gè)點(diǎn).可設(shè)一個(gè)函數(shù)Function distPtoP計(jì)算點(diǎn)與點(diǎn)之間的距離為
再對(duì)圓弧的角度編程,并區(qū)分出不同的象限的各種情況,即其中:t2為終止角;t1為起始角;S為圓弧長(zhǎng).圓弧坐標(biāo)為:x=R·cos t1;y=R·sin t1.
調(diào)用Function distPtoP和Function Get Angle計(jì)算半徑和起始角、終止角.利用VB中的Circle命令畫(huà)圓弧
網(wǎng)格劃分的一般思路是將四邊形的四條邊獨(dú)立看待,可以用插入等比數(shù)列等方式控制網(wǎng)格的疏密.
平面有限元網(wǎng)格劃分時(shí),可按圖1所示的流程進(jìn)行.
圖1 有限元網(wǎng)格劃分的流程圖Fig.1 A flow chart for finite element grid dividing
例1 對(duì)于一個(gè)板中開(kāi)了2個(gè)圓孔的構(gòu)件可以劃分為16個(gè)帶弧邊的四邊形和16個(gè)直線(xiàn)四邊形,由于這些子區(qū)域有的完全相同,可使用鏡像功能.靠近圓孔處的單元要密集些,然后依次向外部擴(kuò)大,輸入一個(gè)比例系數(shù)為1.1的等比數(shù)列即可,再輸入劃分的單元數(shù),以1/4圓所在的弧邊四邊形為例,在兩個(gè)方向上分別為6和6,單擊劃分就可得到36個(gè)按規(guī)則排列的子單元,如圖2所示.
圖2 含2個(gè)圓孔的平面有限元網(wǎng)格劃分Fig.2 Finite element grid for a plane with 2 holes
例2 如圖3所示的曲柄可視為由若干個(gè)雙邊為弧邊的四邊形和4個(gè)一邊為弧邊的四邊形組合.與上例相同,輸入所需的比例系數(shù)和劃分單元數(shù)即可劃分得到了曲柄的有限元網(wǎng)格.
圖3 曲柄的平面有限元網(wǎng)格劃分Fig.3 Finite element grid for a crank
參照當(dāng)前比較成熟的四邊形網(wǎng)格劃分算法,針對(duì)邊界較復(fù)雜的小型構(gòu)件并采用Visual Basic程序語(yǔ)言開(kāi)發(fā)出了四邊形單元的有限元網(wǎng)格劃分程序,可以對(duì)復(fù)雜的平面區(qū)域進(jìn)行劃分,得到較高質(zhì)量的四邊形單元.生成的結(jié)構(gòu)數(shù)據(jù)合理,為下一步進(jìn)行平面有限元的計(jì)算奠定了基礎(chǔ).
[1]吳淑芳,張樹(shù)仁.四叉樹(shù)法生成網(wǎng)格的幾何復(fù)雜度[J].兵工學(xué)報(bào),2000,21(3):282-285.
[2]關(guān)振群,宋超,顧元憲,等.有限元網(wǎng)格生成方法研究的新進(jìn)展[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2003,15(1):1-11.
[3]杜平安.有限元網(wǎng)格劃分的基本原則[J].機(jī)械設(shè)計(jì)與制造,2000(1):34-36.
[4]馬新武,趙國(guó)群.全自動(dòng)四邊形有限元網(wǎng)格生成—I:區(qū)域分解法[J].塑性工程學(xué)報(bào),2007,14(4):105-109.
[5]趙熠,趙建軍.前沿法生成四邊形網(wǎng)格的改進(jìn)方法[J].機(jī)械工程學(xué)報(bào),2007,18(3):308-312.
[6]BLACKER T D,STEPHENSON M B.Paving:a new approach to automated quadrilateral mesh generation[J].Int J Numer Mesh Engng,1991,32:811-847.
[7]SARRATE J,HUERTA A.Efficient unstructured quadrilateral mesh generation[J].Int J Number Mesh Engng,2000,49(9):1327-1350.
[8]蘇金明.用Visual Basic開(kāi)發(fā)交互式CAD系統(tǒng)[M].北京:電子工業(yè)出版社,2003.
Study of a practical dividing method for the plane finite element grid
PU Jun-ping,YAO Yu-long
(College of Civil Engineering & Architecture,Zhejiang University of Technology,Hangzhou 310032,China)
The generation of grid is the precondition and foundation of plane finite element calculation,but poor units will be generated when there are irregular boundaries,which will affects the calculation results.Referring matured quadrilateral grid generation methods,a Visual basic program based on structured finite element grid portioning is developed,which will generate regular quadrilateral grid.It has the characteristics such as:easy to use,grid density and united quantity could be controlled,the data generated are brief and practicable,especially for finite element analysis of small structures.
grid generation;quadrilateral element;grid density;Visual Basic language
TP391
A
1006-4303(2012)01-0088-04
2010-10-27
蒲軍平(1962—),男,新疆烏魯木齊人,教授,博士后,主要從事結(jié)構(gòu)工程、結(jié)構(gòu)動(dòng)力學(xué)研究,E-mail:pjp@zjut.edu.cn.
(
劉 巖)