孫洪鐵
(中國(guó)天辰工程有限公司,天津 300400)
有限單元法是一種獲得工程問(wèn)題近似解的數(shù)值分析方法,它通過(guò)離散的有限個(gè)單元集合體來(lái)代替真實(shí)的結(jié)構(gòu),通過(guò)對(duì)單元的分析,進(jìn)而得到整個(gè)結(jié)構(gòu)的近似解,在工程設(shè)計(jì)中,這種近似解已能滿足工程的需要,有足夠的精確度。有限單元法是最初結(jié)構(gòu)力學(xué)中桿系結(jié)構(gòu)的矩陣位移法的一種推廣,并應(yīng)用于二維和三維問(wèn)題。然而,與桿系結(jié)構(gòu)不同,二維和三維實(shí)體沒(méi)有明顯的聯(lián)結(jié)點(diǎn),這就需要建立許多人為的節(jié)點(diǎn),將連續(xù)體離散化為許多任意形狀的單元,用此方法,連續(xù)體便可用有限個(gè)自由度的系統(tǒng)來(lái)近似表示,并求得其近似解。
從求解未知量角度來(lái)看,有限元法可分為應(yīng)力場(chǎng),位移場(chǎng)及混合場(chǎng)三種場(chǎng)變量,本文以應(yīng)用較為普遍的位移場(chǎng)為例來(lái)闡述有限元的分析過(guò)程,其主要分析步驟為:連續(xù)體的離散化,選擇位移函數(shù),有限元模型的建立,集合離散化單元形成系統(tǒng)方程組,求解節(jié)點(diǎn)位移及通過(guò)位移計(jì)算單元的應(yīng)變和應(yīng)力。
離散化主要目標(biāo)是將物體分成充分小的單元,使得簡(jiǎn)單的位移模型即能足夠近似的表示精確解,在這個(gè)過(guò)程中,將人為通過(guò)網(wǎng)格劃分出有限個(gè)單元,單元與單元之間以節(jié)點(diǎn)相連。以二維問(wèn)題為例,節(jié)點(diǎn)的設(shè)置要遵循以下幾點(diǎn)要求:集中載荷處,分布載荷突變處,幾何形狀不連續(xù)處,材料性質(zhì)突變處,幾何形狀的凹角處等,單元的形狀主要有三角形,矩形等,三角形要求不要出現(xiàn)鈍角;矩形盡量不要采用長(zhǎng)寬比較大的矩形,否則都會(huì)影響近似解的精度。
有限元法的基本原理是分塊近似,即把所要研究的區(qū)域,分割成有限個(gè)子區(qū)域,然后假設(shè)一些比較簡(jiǎn)單的函數(shù)來(lái)表示每個(gè)子區(qū)域中的解,這些假定的函數(shù)被稱為位移函數(shù),位移場(chǎng)或位移模式。
一般多項(xiàng)式被作為位移函數(shù),如二維位移函數(shù)一般為:
不僅因?yàn)槠湓谖⒎峙c積分計(jì)算中比較容易,而且任一階的多項(xiàng)式可以近似表示真實(shí)解,我們知道無(wú)限多次多項(xiàng)式與精確解相對(duì)應(yīng),我們?cè)诓煌碾A次將其截?cái)?,就得到不同程度的近似?如圖1a)~圖1c)所示);其次有限元法作為一種數(shù)值計(jì)算,所選擇的位移函數(shù)一定要保證其收斂或趨向于問(wèn)題的精確解,必須要滿足以下三個(gè)條件:1)位移函數(shù)在單元內(nèi)必須連續(xù),并且相鄰單元的位移必須協(xié)調(diào)。單元內(nèi)的連續(xù)性,通過(guò)選擇具備連續(xù)性的多項(xiàng)式即能滿足,相鄰單元之間的位移協(xié)調(diào),要求單元之間不能開(kāi)裂、重疊,這要求單元邊界的位移僅與該邊界上各節(jié)點(diǎn)的位移有關(guān)。2)位移函數(shù)必須包含單元的剛體位移,即單元一定要產(chǎn)生但不引起應(yīng)力的那部分位移,如平面應(yīng)力單元會(huì)在平面內(nèi)任意方向產(chǎn)生均勻移動(dòng)和轉(zhuǎn)動(dòng)而不引起應(yīng)變。在多項(xiàng)式位移函數(shù)中,常數(shù)項(xiàng)提供單元的剛體位移。3)位移函數(shù)必須包含單元的常應(yīng)變狀態(tài),從物理上我們可以理解常應(yīng)變狀態(tài)的必要性,設(shè)想表示結(jié)構(gòu)集合體的單元越來(lái)越多,則在極限情況下,每個(gè)單元趨近于一個(gè)非常小的尺寸,單元內(nèi)的應(yīng)變接近為常量,所以位移函數(shù)中必須包含常應(yīng)變,才能使計(jì)算收斂于正確解,多項(xiàng)式一次項(xiàng)提供了單元的常應(yīng)變。以上1)條稱為有限單元的協(xié)調(diào)性;2),3)條稱做有限單元的完備性,如單元既完備又協(xié)調(diào),則收斂是單調(diào)的,即以某種模來(lái)度量的分析結(jié)果的精度會(huì)隨著單元數(shù)目的增加而不斷提高,如果只具備完備但不協(xié)調(diào),則分析結(jié)果在極限時(shí)也可能收斂于“精確”結(jié)果,但一般收斂是不單調(diào)的,但在實(shí)際工程中,非協(xié)調(diào)元有時(shí)比與它密切相關(guān)的協(xié)調(diào)單元要好,原因在于采用的近似解的性質(zhì),我們都知道,利用假定的位移函數(shù)得到的近似結(jié)構(gòu)比實(shí)際結(jié)構(gòu)要更剛一些,但由于允許單元分離,重疊使這種近似結(jié)構(gòu)變?nèi)崃耍@兩種影響相互抵消,常常得到好的結(jié)果,比如非協(xié)調(diào)元在板單元中的運(yùn)用。選擇位移函數(shù)的另一個(gè)因素,即該函數(shù)與局部坐標(biāo)系的方位無(wú)關(guān),即在任何一組相對(duì)于單元來(lái)說(shuō)的方向固定的荷載作用下,單元的反應(yīng)(指在和單元一起移動(dòng)的坐標(biāo)系中的單元應(yīng)變能或應(yīng)變)不依賴于單元本身及它的荷載在全局坐標(biāo)xy中的方向,單元的這種性質(zhì)叫做單元的幾何各向同性或幾何不變性,對(duì)于一般的多項(xiàng)式位移函數(shù),按帕斯卡三角形(如圖1d)所示)對(duì)稱選取即可。
圖1 位移函數(shù)
作為位移場(chǎng)中有限單元的建立,我們所要求的是連接有限個(gè)單元體的節(jié)點(diǎn)的位移與節(jié)點(diǎn)力之間的轉(zhuǎn)化關(guān)系,即確定單元的剛度矩陣。有限元矩陣的建立一般分廣義坐標(biāo)有限模型和等參有限元模型,前者在2.2節(jié)中位移函數(shù)式(1)我們已經(jīng)看到,其中α1,α2,α3,…,αm被稱為廣義坐標(biāo),通過(guò)廣義坐標(biāo)將單元內(nèi)任一點(diǎn)的位移與單元節(jié)點(diǎn)位移相聯(lián)系,并且通過(guò)一個(gè)非奇異矩陣來(lái)求得以節(jié)點(diǎn)位移及節(jié)點(diǎn)坐標(biāo)表示的單元內(nèi)任意一點(diǎn)的位移。而等參有限元模型建立的基本思想是:直接通過(guò)使用插值函數(shù)(或稱形函數(shù))來(lái)得到單元內(nèi)任意一點(diǎn)的位移和單元節(jié)點(diǎn)位移之間的關(guān)系。在實(shí)際分析中,使用等參有限元有時(shí)更為有效,可以避免在廣義坐標(biāo)計(jì)算中出現(xiàn)奇異矩陣的可能性,也減少了矩陣運(yùn)算的次數(shù);同時(shí)二次或高階等參單元既可以模擬直邊也可以模擬曲邊,在模擬曲線邊界的結(jié)構(gòu)時(shí)非常有用。然而廣義坐標(biāo)有限元模型的建立,能更好的讓我們加深對(duì)有限元法的理解,因此本文以廣義坐標(biāo)法為例來(lái)介紹有限單元模型的建立,即求解單元?jiǎng)偠染仃嚨倪^(guò)程。從推導(dǎo)方法來(lái)看,分為三類:1)直接法,該方法易于理解,適用于簡(jiǎn)單的問(wèn)題。2)變分法,把有限元?dú)w結(jié)為求泛函的極值問(wèn)題,比如固體力學(xué)中的最小勢(shì)能原理的應(yīng)用,它使有限元法建立在更加堅(jiān)實(shí)的數(shù)學(xué)基礎(chǔ)上,擴(kuò)大了有限元法的應(yīng)用范圍。3)加權(quán)余數(shù)法,不需要利用泛函的概念,而直接從基本的微分方程出發(fā),求出近似解,對(duì)于不存在泛函的工程領(lǐng)域,提供了有效的解決方法。
以下我們將通過(guò)比較直觀易于理解的直接法來(lái)了解一下有限元模型的建立過(guò)程,以彈性力學(xué)平面應(yīng)力問(wèn)題為例(見(jiàn)圖2):首先假定線性多項(xiàng)式為位移函數(shù),使廣義坐標(biāo)α個(gè)數(shù)等于單元節(jié)點(diǎn)位移個(gè)數(shù):
可以寫成:
可簡(jiǎn)寫為:
此為單元內(nèi)部點(diǎn)位移與廣義坐標(biāo)之間的轉(zhuǎn)化關(guān)系,為了求出節(jié)點(diǎn)位移與單元內(nèi)部點(diǎn)位移的關(guān)系,需要求出節(jié)點(diǎn)位移{δ}(e)與{α}的轉(zhuǎn)換關(guān)系。
三角形單元的節(jié)點(diǎn)位移可以表示為:
{δ}(e)= [{δ1},{δ2},{δ3}]T=[υ1v1υ2v2υ3v3]T。
三角形單元的節(jié)點(diǎn)的力向量可以表示為:
{F}(e)= [{F1},{F2},{F3}]T=[u1υ1u2υ2u3υ3]T。
我們要找到{δ}(e)與{F}(e)的轉(zhuǎn)換關(guān)系,即{F}(e)=[k](e)×{δ}(e),其中,[k](e)為單元的剛度矩陣。以節(jié)點(diǎn)的水平分量為例代入式(2),可以寫成:
節(jié)點(diǎn)1:υ1=α1+α2x1+α3y1;節(jié)點(diǎn)2:υ2=α1+α2x2+α3y2;節(jié)點(diǎn)3:υ3=α1+α2x3+α3y3。
由上面三個(gè)方程,易得 α1,α2,α3;同理,通過(guò)節(jié)點(diǎn)豎直分力,可得α4,α5,α6,最終可以得到一個(gè)廣義坐標(biāo)與節(jié)點(diǎn)位移的關(guān)系式:
其中,[A]為一個(gè)只與三角形單元節(jié)點(diǎn)坐標(biāo)有關(guān)的矩陣。將式(5)代入式(4),得:
至此單元內(nèi)部任意點(diǎn)位移已由節(jié)點(diǎn)坐標(biāo)及節(jié)點(diǎn)位移表示。接下來(lái)可以求解單元的應(yīng)變與位移的關(guān)系,在彈性力學(xué)平面問(wèn)題中由:
求得:
轉(zhuǎn)換矩陣[B]是一個(gè)僅僅與三角形單元的幾何性質(zhì)有關(guān)的常量,被稱為幾何矩陣,由上式我們確定單元位移場(chǎng)的應(yīng)變狀態(tài)。
為引入單元材料性質(zhì)的影響,需要通過(guò)應(yīng)力與應(yīng)變的本構(gòu)方程,求解單元應(yīng)力:
轉(zhuǎn)換矩陣[D]被稱為彈性矩陣,對(duì)平面應(yīng)力問(wèn)題:
接下來(lái),由單元應(yīng)力推算節(jié)點(diǎn)力,需要用到平衡方程,在有限元法中通常用虛功原理來(lái)代替平衡方程。在平面應(yīng)力問(wèn)題中虛功原理表達(dá)式為:
如圖3所示中,令實(shí)際受力狀態(tài)在虛設(shè)位移狀態(tài)上作虛功,由{ε*}=[B]{δ*}(e)得:
代入式(10)中,由于{δ*}為任意的,整理后得:
因在常應(yīng)變?nèi)切螁卧校跙],{σ}均為常量,則:
結(jié)合式(7)及式(8)得:
其中,[k](e)為單元的剛度矩陣。至此,我們已經(jīng)完成了單元分析,為單元的集成提供了必要的條件。
圖2 平面應(yīng)力三角形單元節(jié)點(diǎn)位移及節(jié)點(diǎn)力
圖3 平面應(yīng)力三角形單元節(jié)點(diǎn)力虛設(shè)位移
單元?jiǎng)偠燃傻姆椒ɑ诟饔邢迒卧w連接節(jié)點(diǎn)處的協(xié)調(diào)性要求,即有限單元在相互連接的節(jié)點(diǎn)處,必須具有相同的位移,這里的位移為廣義位移,它可以是位移、轉(zhuǎn)角。出于這種考慮,我們采用直接剛度法建立整體剛度矩陣[K],因?yàn)閱卧獎(jiǎng)偠染仃嚕踜](e)與整體剛度矩陣[K]的階數(shù)并不相同,需要將單元?jiǎng)偠染仃嚁U(kuò)大成與整體剛度矩陣同階數(shù),擴(kuò)大后的矩陣叫做單元的貢獻(xiàn)矩陣,它們表示每個(gè)單元單獨(dú)變形時(shí)對(duì)整體剛度矩陣提供的貢獻(xiàn)。然后根據(jù)局部編碼和整體編碼的對(duì)應(yīng)關(guān)系,進(jìn)行單元?jiǎng)偠染仃嚨寞B加,節(jié)點(diǎn)荷載也采用同樣的方式進(jìn)行疊加。在有限元法中,通常在整體剛度矩陣[K]形成后,引入已知的邊界條件,而后形成可求解的有限元系統(tǒng)方程組:{F}=[K]{δ},其中,[K]為整體剛度矩陣;{F}為結(jié)構(gòu)的節(jié)點(diǎn)力;{δ}為結(jié)構(gòu)節(jié)點(diǎn)位移。
有限元系統(tǒng)方程組一旦建立,我們可以根據(jù)適用于計(jì)算機(jī)處理的方程組的解法,如高斯消元法等,求解出結(jié)構(gòu)各有限單元節(jié)點(diǎn)處的位移,然后按式(7)及式(8)即可算出有限單元任意點(diǎn)處應(yīng)力及應(yīng)變分量,然而我們必須指出一點(diǎn),單元的應(yīng)力并不能保證單個(gè)單元的平衡條件,即在單元的節(jié)點(diǎn)處內(nèi)應(yīng)力和所加載荷是不一定平衡的,我們通常用應(yīng)力平均值來(lái)代表單元的應(yīng)力,最常用的平均法是采用單元形心處的應(yīng)力。
有限單元法現(xiàn)今已成為工程師進(jìn)行結(jié)構(gòu)分析的有效工具,并可廣泛的應(yīng)用于各個(gè)領(lǐng)域,然而有限單元的分析過(guò)程,基本都是通過(guò)計(jì)算機(jī)編程來(lái)完成的,本文力求通過(guò)對(duì)有限單元的概述,使工程設(shè)計(jì)人員對(duì)其基本概念有個(gè)宏觀的認(rèn)識(shí),以便使工程設(shè)計(jì)人員更好的使用有限元分析軟件來(lái)解決實(shí)際問(wèn)題。
[1]龍馭球.有限元法概論[M].北京:人民教育出版社,1979:381.
[2]R.D.庫(kù)克.有限元分析的概念和應(yīng)用[M].何 窮,程耿東,譯.北京:科技出版社,1981.
[3]K.J.巴特,E.L.威爾遜.有限元分析中的數(shù)值方法[M].林公豫,羅 恩,譯.北京:科技出版社,1985.