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

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

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

可以寫成:
可簡寫為:

此為單元內部點位移與廣義坐標之間的轉化關系,為了求出節點位移與單元內部點位移的關系,需要求出節點位移{δ}(e)與{α}的轉換關系。
三角形單元的節點位移可以表示為:
{δ}(e)= [{δ1},{δ2},{δ3}]T=[υ1v1υ2v2υ3v3]T。
三角形單元的節點的力向量可以表示為:
{F}(e)= [{F1},{F2},{F3}]T=[u1υ1u2υ2u3υ3]T。
我們要找到{δ}(e)與{F}(e)的轉換關系,即{F}(e)=[k](e)×{δ}(e),其中,[k](e)為單元的剛度矩陣。以節點的水平分量為例代入式(2),可以寫成:
節點1:υ1=α1+α2x1+α3y1;節點2:υ2=α1+α2x2+α3y2;節點3:υ3=α1+α2x3+α3y3。
由上面三個方程,易得 α1,α2,α3;同理,通過節點豎直分力,可得α4,α5,α6,最終可以得到一個廣義坐標與節點位移的關系式:

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

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

求得:

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

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

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

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

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

因在常應變三角形單元中[B],{σ}均為常量,則:

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

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

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

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