吳 敏
(西南交通大學,四川 成 都 610000)
有限元分析是求解復雜微分方程的一種非常有效的數值方法[1],在科學研究和工程分析中廣泛運用。它結合計算機技術,具有計算速度快、計算精度高的特點,廣泛運用于力學領域、熱學領域、電磁學領域以及學科交叉領域中。有限元思想起源于1943年Courant[2]解決St.Venant扭轉問題的工作中,Clough[3]在1960年首次提出“有限單元法”的名稱。有限元的核心是將一個連續的物體通過離散化的形式來表示,將連續體離散為若干個單元的組合,每個單元通過節點的形式連接。通過離散化的方式將具有無限自由度的物體用有限自由度的離散體近似表示,從而得到相關的數值解。運用有限元對結構的受力情況進行數值仿真時,以節點位移為基本未知量,首先用節點位移來表示離散體單元內部的應力和應變狀態,然后通過虛功原理構造單元剛度矩陣,接著將單元剛度矩陣組裝成為總體剛度矩陣,并施加邊界條件求出每個節點上的位移值,最后將節點上的位移值帶入到每個單元中,求解出每個單元的應力和應變狀態。
目前關于有限元的文獻中,外加載荷的形式普遍采用的是力載荷,幾乎沒有關于外載荷形式為位移載荷的有限元分析原理的文獻,針對這種現象,本文分別推導了外加載荷為位移載荷的形式下有限元原理,并且使用MATLAB 軟件分別編寫出了外載荷為力載荷和位移載荷的有限元程序。對比兩種外載荷形式下的數值仿真結果發現:外載荷為位移載荷計算出的仿真結果與外載荷為力載荷計算出的仿真結果基本一致。由此可以拓展有限元的分析思路:當外加力載荷不容易確定時,通過測量位移載荷的大小來對結構進行有限元分析,同樣可以得到結構的有限元受力結果。

圖1 結構的有限元受力分析步驟
在對結構的受力情況進行有限元分析時,以節點位移為基本未知量,主要分析步驟采用的是“總-分-總-分”的形式[4],即首先將結構離散化為若干個單元,然后再求出單元剛度矩陣,接著將若干個單元剛度矩陣組裝為一個整體剛度矩陣,并且施加邊界約束和外加載荷,求出節點位移,最后將求出的節點位移帶入到每一個單元中,求解出每一個單元的應變和應力情況。有限元受力分析的流程圖如圖1 所示。在本文中,以一個二維平面的連續體為例,來推導有限元的分析步驟[5]。其中,單元類型采用的是的目前廣泛使用的三節點三角形單元,如圖2 所示。

圖2 三節點三角形單元
將一個二維的連續體離散為具有若干個三節點三角形單元,通過離散化處理,將具有無限自由度的連續體近似等效為有限自由度的離散體。圖3(a)和(b)分別表示一個二維的連續體和三節點三角形單元表示的離散體。

圖3 二維平面的連續體和離散體對應的結構示意圖
將結構經過離散化處理后,結構由若干個離散的單元表示。在單元分析的過程中,將用單元的節點位移分別表示單元內部位移、單元應變和單元應力,并且求出單元剛度矩陣。
2.2.1 單元內部位移
在每一個單元中,引入形函數,則單元內部任意一點的位移可以通過單元上的節點位移來表示。

式(1)中,ue=[uiviujvjukvk]T表示單元的節點位移,a=[uv]T表示單元內部任意一點的位移,N 表示形函數矩陣,表達形式為:

式(2)中,Ni、Nj、Nk是與節點坐標和位移待求點坐標有關的函數,稱為形函數。在這里,形函數采用一次多項式的形式,表達形式為:

式(3)中,D 代表三角形單元的面積,x、y 表示單元任意一點的坐標,其他系數的表示形式為:

式(4)中,“(i,j,k)”的意思是依次將下標做輪換:i→j,j→k,k→i。這樣,所有的待定參數都可以通過節點位置來表示。
2.2.2 單元應變
在確定了單元內任意一點的位移之后,通過幾何方程即可確定單元內的應變:

式(5)中,∈為單元應變矩陣,B 為形函數導數矩。在三節點三角形單元中,B 的表示形式為:

2.2.3 單元應力
求出單元應變以后,利用物理方程可以求出單元內的應力分布:

在二維平面情況下,結構可以分為平面應力問題和平面應變問題。對于平面應力問題,D 可以表示為:

對于平面應變問題,D 可以表示為:

式(8)和式(9)中,E 為彈性模量,μ 為剪切模量。
2.2.4 求解單元剛度矩陣
在求出應變和應力的表達方式以后,通過虛位移原理可以得到單元剛度矩陣和單元節點力列陣。虛位移原理可以表示為:

式(10)中,f 和F 分別表示體積力和面積力。在離散化結構中,單元內的虛位移原理可以表示為:

經過化簡后,可以得到單元內部的有限元方程:

在整體分析的過程中,首先將單元剛度矩陣組裝成整體剛度矩陣,然后施加邊界條件和外加載荷,最后計算求解出節點位移。
2.3.1 求解整體剛度矩陣
在求出單元剛度矩陣和單元節點力列陣后,將若干個單元的單元剛度矩陣和節點力列陣進行組裝,即:

式(13)就是結構受力分析的有限元方程,將其簡化書寫為:

式(14)中,K 為整體剛度矩陣,u 為整體位移列陣,P為整體節點力列陣。將有限元方程展開為矩陣形式,可以寫為:

2.3.2 施加邊界條件和外載荷
邊界條件通常是限制結構在某些位置上的位移,表達形式為:

外載荷的加載形式有位移加載和力加載兩種形式。兩種施加方式的求解過程有所不同,下面分別推導外載荷形式為力載荷和位移載荷下的有限元方程的求解過程。
(1)外載荷形式為力加載的有限元求解過程
當外載荷形式為力載荷時,力載荷的加載形式可以表示為:

將邊界條件和力載荷加入到整體分析的有限元方程中并用矩陣形式表示,可以得到:

式(18)中,u1和 P2是已知的參數,u2和 u3為待求的節點位移值。此時,可以將有限元方程簡化為:

由式(19)即可求得待求的節點位移u2和u3。
(2)外載荷形式為位移載荷的有限元求解過程
當外載荷形式為位移載荷時,位移載荷的加載形式可以表示為:

將邊界條件和位移載荷加入到整體分析的有限元方程中并用矩陣形式表示,可以得到:

式(21)中,u1和 u2是已知的參數,u3和 P2分別為待求的節點位移值和節點載荷值。此時,可以將有限元方程簡化為:

由式(22)即可求得待求的節點位移u3,然后將求出的節點位移值帶入到式(21)中即可求得力載荷P2。
在求出離散體的節點位移以后,首先通過式(5)求解出每個單元的應變分布,然后通過式(7)求解出每個單元的應力分布。
在數值算例中,我們運用MATLAB 軟件編寫出力載荷和位移載荷下的有限元程序,在MATLAB 軟件中完成了力載荷和位移載荷下的有限元仿真,并且得到相關的數值結果。
如圖4 所示為二維平面的矩形板。矩形板的長為3m,寬為 1m,彈性模量為 E=1×107Pa,泊松比 μ=0.3,在矩形板右上方的頂點處施加的力載荷大小為F=1×105N。使用三節點三角形單元,并采用平面應力假設,求解矩形板每個節點上的位移。

圖4 力載荷作用下的矩形板及其有限元幾何模型
受力載荷作用下的有限元分析步驟如下所示:
(1)將矩形板進行離散化處理;
(2)使用式(12)求解出每一個單元的單元剛度矩陣;
(3)使用式(13)組裝整體剛度矩陣;
(4)處理邊界條件和力載荷
矩形板的邊界條件為:
u7=v7=u8=v8=0;
外加的節點力載荷為:
Py1=-1×104N;
(5)求解未知的節點位移
運用式(18)和式(19),求出節點位移列陣為:
u=[0.0069-0.0329-0.0058-0.0315 0.0057-0.018-0.0052-0.0174 0.0035-0.0064-0.0034-0.0059 0 0 0 0]T;
(6)求解單元內的應力狀態
將節點位移帶入到每一個單元內,可求出力載荷情況下單元1~6 中的應力狀態,具體數值如表1 所示。
如圖5 所示為一個二維平面的方形板。矩形板的長為3m,寬為 1m,彈性模量為 E=1×107Pa,泊松比 μ=0.3,在矩形板右上方的頂點處施加的位移載荷大小為u=-0.0329m。使用三節點三角形單元,并采用平面應力假設,求解矩形板每個節點上的位移和每個單元的應力值。

圖5 位移載荷作用下的矩形板及其有限元幾何模型
受位移載荷作用下的有限元分析步驟如下所示:
(1)將矩形板進行離散化處理;
(2)使用式(12)求解出每一個單元的單元剛度矩陣;
(3)使用式(13)組裝整體剛度矩陣;

表1 力載荷下各單元的應力值(Pa)

表2 位移載荷下各單元的應力值(Pa)
(4)處理邊界條件和位移載荷
矩形板的邊界條件為:
u7=v7=u8=v8=0;
外加的節點位移載荷為:
v1=-0.0329m;
(5)求解未知的節點位移和外加力載荷
使用式(21)和式(22)求解未知的位移,求出的節點位移列陣為:
u=[0.0069-0.0329-0.0058-0.0315 0.0057-0.018-0.0052-0.0174 0.0035-0.0064-0.0034-0.0059 0 0 0 0]T;
使用式(21),求得施加位移載荷時等效的外加載荷,其值為:
Py1=-0.993×103N;
(6)求解各單元的應力情況
將節點位移帶入到每一個單元內,可求出位移載荷情況下單元1~6 中的應力狀態,具體數值如表2 所示。
使用力載荷加載計算出的位移結果和使用位移載荷計算出的位移結果幾乎一致,兩種外加載荷計算出來的應力結果也大致相同。說明使用兩種不同的外加載荷進行有限元分析得到的數值結果是一致的,驗證了運用位移載荷進行有限元分析的正確性。
本文分析了外加載荷分別為力載荷和位移載荷情況下的有限元原理,并且使用MATLAB 軟件編寫了兩種外加載荷在二維平面應力問題下的有限元程序,驗證了運用力載荷和位移載荷對結構進行有限元分析時仿真結果的一致性。為結構的有限元分析擴展了思路。當結構的外加力載荷不容易得到時,通過測量位移載荷的大小也可以對結構進行有限元分析。