孫洪鐵
(中國天辰工程有限公司,天津 300400)
有限單元法是一種獲得工程問題近似解的數(shù)值分析方法,它通過離散的有限個單元集合體來代替真實的結(jié)構(gòu),通過對單元的分析,進而得到整個結(jié)構(gòu)的近似解,在工程設(shè)計中,這種近似解已能滿足工程的需要,有足夠的精確度。有限單元法是最初結(jié)構(gòu)力學中桿系結(jié)構(gòu)的矩陣位移法的一種推廣,并應用于二維和三維問題。然而,與桿系結(jié)構(gòu)不同,二維和三維實體沒有明顯的聯(lián)結(jié)點,這就需要建立許多人為的節(jié)點,將連續(xù)體離散化為許多任意形狀的單元,用此方法,連續(xù)體便可用有限個自由度的系統(tǒng)來近似表示,并求得其近似解。
從求解未知量角度來看,有限元法可分為應力場,位移場及混合場三種場變量,本文以應用較為普遍的位移場為例來闡述有限元的分析過程,其主要分析步驟為:連續(xù)體的離散化,選擇位移函數(shù),有限元模型的建立,集合離散化單元形成系統(tǒng)方程組,求解節(jié)點位移及通過位移計算單元的應變和應力。
離散化主要目標是將物體分成充分小的單元,使得簡單的位移模型即能足夠近似的表示精確解,在這個過程中,將人為通過網(wǎng)格劃分出有限個單元,單元與單元之間以節(jié)點相連。以二維問題為例,節(jié)點的設(shè)置要遵循以下幾點要求:集中載荷處,分布載荷突變處,幾何形狀不連續(xù)處,材料性質(zhì)突變處,幾何形狀的凹角處等,單元的形狀主要有三角形,矩形等,三角形要求不要出現(xiàn)鈍角;矩形盡量不要采用長寬比較大的矩形,否則都會影響近似解的精度。
有限元法的基本原理是分塊近似,即把所要研究的區(qū)域,分割成有限個子區(qū)域,然后假設(shè)一些比較簡單的函數(shù)來表示每個子區(qū)域中的解,這些假定的函數(shù)被稱為位移函數(shù),位移場或位移模式。
一般多項式被作為位移函數(shù),如二維位移函數(shù)一般為:

不僅因為其在微分與積分計算中比較容易,而且任一階的多項式可以近似表示真實解,我們知道無限多次多項式與精確解相對應,我們在不同的階次將其截斷,就得到不同程度的近似解(如圖1a)~圖1c)所示);其次有限元法作為一種數(shù)值計算,所選擇的位移函數(shù)一定要保證其收斂或趨向于問題的精確解,必須要滿足以下三個條件:1)位移函數(shù)在單元內(nèi)必須連續(xù),并且相鄰單元的位移必須協(xié)調(diào)。單元內(nèi)的連續(xù)性,通過選擇具備連續(xù)性的多項式即能滿足,相鄰單元之間的位移協(xié)調(diào),要求單元之間不能開裂、重疊,這要求單元邊界的位移僅與該邊界上各節(jié)點的位移有關(guān)。2)位移函數(shù)必須包含單元的剛體位移,即單元一定要產(chǎn)生但不引起應力的那部分位移,如平面應力單元會在平面內(nèi)任意方向產(chǎn)生均勻移動和轉(zhuǎn)動而不引起應變。在多項式位移函數(shù)中,常數(shù)項提供單元的剛體位移。3)位移函數(shù)必須包含單元的常應變狀態(tài),從物理上我們可以理解常應變狀態(tài)的必要性,設(shè)想表示結(jié)構(gòu)集合體的單元越來越多,則在極限情況下,每個單元趨近于一個非常小的尺寸,單元內(nèi)的應變接近為常量,所以位移函數(shù)中必須包含常應變,才能使計算收斂于正確解,多項式一次項提供了單元的常應變。以上1)條稱為有限單元的協(xié)調(diào)性;2),3)條稱做有限單元的完備性,如單元既完備又協(xié)調(diào),則收斂是單調(diào)的,即以某種模來度量的分析結(jié)果的精度會隨著單元數(shù)目的增加而不斷提高,如果只具備完備但不協(xié)調(diào),則分析結(jié)果在極限時也可能收斂于“精確”結(jié)果,但一般收斂是不單調(diào)的,但在實際工程中,非協(xié)調(diào)元有時比與它密切相關(guān)的協(xié)調(diào)單元要好,原因在于采用的近似解的性質(zhì),我們都知道,利用假定的位移函數(shù)得到的近似結(jié)構(gòu)比實際結(jié)構(gòu)要更剛一些,但由于允許單元分離,重疊使這種近似結(jié)構(gòu)變?nèi)崃耍@兩種影響相互抵消,常常得到好的結(jié)果,比如非協(xié)調(diào)元在板單元中的運用。選擇位移函數(shù)的另一個因素,即該函數(shù)與局部坐標系的方位無關(guān),即在任何一組相對于單元來說的方向固定的荷載作用下,單元的反應(指在和單元一起移動的坐標系中的單元應變能或應變)不依賴于單元本身及它的荷載在全局坐標xy中的方向,單元的這種性質(zhì)叫做單元的幾何各向同性或幾何不變性,對于一般的多項式位移函數(shù),按帕斯卡三角形(如圖1d)所示)對稱選取即可。

圖1 位移函數(shù)
作為位移場中有限單元的建立,我們所要求的是連接有限個單元體的節(jié)點的位移與節(jié)點力之間的轉(zhuǎn)化關(guān)系,即確定單元的剛度矩陣。有限元矩陣的建立一般分廣義坐標有限模型和等參有限元模型,前者在2.2節(jié)中位移函數(shù)式(1)我們已經(jīng)看到,其中α1,α2,α3,…,αm被稱為廣義坐標,通過廣義坐標將單元內(nèi)任一點的位移與單元節(jié)點位移相聯(lián)系,并且通過一個非奇異矩陣來求得以節(jié)點位移及節(jié)點坐標表示的單元內(nèi)任意一點的位移。而等參有限元模型建立的基本思想是:直接通過使用插值函數(shù)(或稱形函數(shù))來得到單元內(nèi)任意一點的位移和單元節(jié)點位移之間的關(guān)系。在實際分析中,使用等參有限元有時更為有效,可以避免在廣義坐標計算中出現(xiàn)奇異矩陣的可能性,也減少了矩陣運算的次數(shù);同時二次或高階等參單元既可以模擬直邊也可以模擬曲邊,在模擬曲線邊界的結(jié)構(gòu)時非常有用。然而廣義坐標有限元模型的建立,能更好的讓我們加深對有限元法的理解,因此本文以廣義坐標法為例來介紹有限單元模型的建立,即求解單元剛度矩陣的過程。從推導方法來看,分為三類:1)直接法,該方法易于理解,適用于簡單的問題。2)變分法,把有限元歸結(jié)為求泛函的極值問題,比如固體力學中的最小勢能原理的應用,它使有限元法建立在更加堅實的數(shù)學基礎(chǔ)上,擴大了有限元法的應用范圍。3)加權(quán)余數(shù)法,不需要利用泛函的概念,而直接從基本的微分方程出發(fā),求出近似解,對于不存在泛函的工程領(lǐng)域,提供了有效的解決方法。
以下我們將通過比較直觀易于理解的直接法來了解一下有限元模型的建立過程,以彈性力學平面應力問題為例(見圖2):首先假定線性多項式為位移函數(shù),使廣義坐標α個數(shù)等于單元節(jié)點位移個數(shù):

可以寫成:
可簡寫為:

此為單元內(nèi)部點位移與廣義坐標之間的轉(zhuǎn)化關(guān)系,為了求出節(jié)點位移與單元內(nèi)部點位移的關(guān)系,需要求出節(jié)點位移{δ}(e)與{α}的轉(zhuǎn)換關(guān)系。
三角形單元的節(jié)點位移可以表示為:
{δ}(e)= [{δ1},{δ2},{δ3}]T=[υ1v1υ2v2υ3v3]T。
三角形單元的節(jié)點的力向量可以表示為:
{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é)點的水平分量為例代入式(2),可以寫成:
節(jié)點1:υ1=α1+α2x1+α3y1;節(jié)點2:υ2=α1+α2x2+α3y2;節(jié)點3:υ3=α1+α2x3+α3y3。
由上面三個方程,易得 α1,α2,α3;同理,通過節(jié)點豎直分力,可得α4,α5,α6,最終可以得到一個廣義坐標與節(jié)點位移的關(guān)系式:

其中,[A]為一個只與三角形單元節(jié)點坐標有關(guān)的矩陣。將式(5)代入式(4),得:

至此單元內(nèi)部任意點位移已由節(jié)點坐標及節(jié)點位移表示。接下來可以求解單元的應變與位移的關(guān)系,在彈性力學平面問題中由:

求得:

轉(zhuǎn)換矩陣[B]是一個僅僅與三角形單元的幾何性質(zhì)有關(guān)的常量,被稱為幾何矩陣,由上式我們確定單元位移場的應變狀態(tài)。
為引入單元材料性質(zhì)的影響,需要通過應力與應變的本構(gòu)方程,求解單元應力:

轉(zhuǎn)換矩陣[D]被稱為彈性矩陣,對平面應力問題:

接下來,由單元應力推算節(jié)點力,需要用到平衡方程,在有限元法中通常用虛功原理來代替平衡方程。在平面應力問題中虛功原理表達式為:

如圖3所示中,令實際受力狀態(tài)在虛設(shè)位移狀態(tài)上作虛功,由{ε*}=[B]{δ*}(e)得:

代入式(10)中,由于{δ*}為任意的,整理后得:

因在常應變?nèi)切螁卧校跙],{σ}均為常量,則:

結(jié)合式(7)及式(8)得:

其中,[k](e)為單元的剛度矩陣。至此,我們已經(jīng)完成了單元分析,為單元的集成提供了必要的條件。

圖2 平面應力三角形單元節(jié)點位移及節(jié)點力

圖3 平面應力三角形單元節(jié)點力虛設(shè)位移
單元剛度集成的方法基于各有限單元體連接節(jié)點處的協(xié)調(diào)性要求,即有限單元在相互連接的節(jié)點處,必須具有相同的位移,這里的位移為廣義位移,它可以是位移、轉(zhuǎn)角。出于這種考慮,我們采用直接剛度法建立整體剛度矩陣[K],因為單元剛度矩陣[k](e)與整體剛度矩陣[K]的階數(shù)并不相同,需要將單元剛度矩陣擴大成與整體剛度矩陣同階數(shù),擴大后的矩陣叫做單元的貢獻矩陣,它們表示每個單元單獨變形時對整體剛度矩陣提供的貢獻。然后根據(jù)局部編碼和整體編碼的對應關(guān)系,進行單元剛度矩陣的疊加,節(jié)點荷載也采用同樣的方式進行疊加。在有限元法中,通常在整體剛度矩陣[K]形成后,引入已知的邊界條件,而后形成可求解的有限元系統(tǒng)方程組:{F}=[K]{δ},其中,[K]為整體剛度矩陣;{F}為結(jié)構(gòu)的節(jié)點力;{δ}為結(jié)構(gòu)節(jié)點位移。
有限元系統(tǒng)方程組一旦建立,我們可以根據(jù)適用于計算機處理的方程組的解法,如高斯消元法等,求解出結(jié)構(gòu)各有限單元節(jié)點處的位移,然后按式(7)及式(8)即可算出有限單元任意點處應力及應變分量,然而我們必須指出一點,單元的應力并不能保證單個單元的平衡條件,即在單元的節(jié)點處內(nèi)應力和所加載荷是不一定平衡的,我們通常用應力平均值來代表單元的應力,最常用的平均法是采用單元形心處的應力。
有限單元法現(xiàn)今已成為工程師進行結(jié)構(gòu)分析的有效工具,并可廣泛的應用于各個領(lǐng)域,然而有限單元的分析過程,基本都是通過計算機編程來完成的,本文力求通過對有限單元的概述,使工程設(shè)計人員對其基本概念有個宏觀的認識,以便使工程設(shè)計人員更好的使用有限元分析軟件來解決實際問題。
[1]龍馭球.有限元法概論[M].北京:人民教育出版社,1979:381.
[2]R.D.庫克.有限元分析的概念和應用[M].何 窮,程耿東,譯.北京:科技出版社,1981.
[3]K.J.巴特,E.L.威爾遜.有限元分析中的數(shù)值方法[M].林公豫,羅 恩,譯.北京:科技出版社,1985.