吳 敏

（西南交通大學(xué)，四川 成 都 610000）

1 概述

有限元分析是求解復(fù)雜微分方程的一種非常有效的數(shù)值方法[1]，在科學(xué)研究和工程分析中廣泛運用。它結(jié)合計算機技術(shù)，具有計算速度快、計算精度高的特點，廣泛運用于力學(xué)領(lǐng)域、熱學(xué)領(lǐng)域、電磁學(xué)領(lǐng)域以及學(xué)科交叉領(lǐng)域中。有限元思想起源于1943年Courant[2]解決St.Venant扭轉(zhuǎn)問題的工作中，Clough[3]在1960年首次提出“有限單元法”的名稱。有限元的核心是將一個連續(xù)的物體通過離散化的形式來表示，將連續(xù)體離散為若干個單元的組合，每個單元通過節(jié)點的形式連接。通過離散化的方式將具有無限自由度的物體用有限自由度的離散體近似表示，從而得到相關(guān)的數(shù)值解。運用有限元對結(jié)構(gòu)的受力情況進(jìn)行數(shù)值仿真時，以節(jié)點位移為基本未知量，首先用節(jié)點位移來表示離散體單元內(nèi)部的應(yīng)力和應(yīng)變狀態(tài)，然后通過虛功原理構(gòu)造單元剛度矩陣，接著將單元剛度矩陣組裝成為總體剛度矩陣，并施加邊界條件求出每個節(jié)點上的位移值，最后將節(jié)點上的位移值帶入到每個單元中，求解出每個單元的應(yīng)力和應(yīng)變狀態(tài)。

目前關(guān)于有限元的文獻(xiàn)中，外加載荷的形式普遍采用的是力載荷，幾乎沒有關(guān)于外載荷形式為位移載荷的有限元分析原理的文獻(xiàn)，針對這種現(xiàn)象，本文分別推導(dǎo)了外加載荷為位移載荷的形式下有限元原理，并且使用MATLAB 軟件分別編寫出了外載荷為力載荷和位移載荷的有限元程序。對比兩種外載荷形式下的數(shù)值仿真結(jié)果發(fā)現(xiàn)：外載荷為位移載荷計算出的仿真結(jié)果與外載荷為力載荷計算出的仿真結(jié)果基本一致。由此可以拓展有限元的分析思路：當(dāng)外加力載荷不容易確定時，通過測量位移載荷的大小來對結(jié)構(gòu)進(jìn)行有限元分析，同樣可以得到結(jié)構(gòu)的有限元受力結(jié)果。

圖1 結(jié)構(gòu)的有限元受力分析步驟

2 有限元分析步驟

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

圖2 三節(jié)點三角形單元

2.1 結(jié)構(gòu)離散化

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

圖3 二維平面的連續(xù)體和離散體對應(yīng)的結(jié)構(gòu)示意圖

2.2 單元分析

將結(jié)構(gòu)經(jīng)過離散化處理后，結(jié)構(gòu)由若干個離散的單元表示。在單元分析的過程中，將用單元的節(jié)點位移分別表示單元內(nèi)部位移、單元應(yīng)變和單元應(yīng)力，并且求出單元剛度矩陣。

2.2.1 單元內(nèi)部位移

在每一個單元中，引入形函數(shù)，則單元內(nèi)部任意一點的位移可以通過單元上的節(jié)點位移來表示。

式（1）中，ue=[uiviujvjukvk]T表示單元的節(jié)點位移，a=[uv]T表示單元內(nèi)部任意一點的位移，N 表示形函數(shù)矩陣，表達(dá)形式為：

式（2）中，Ni、Nj、Nk是與節(jié)點坐標(biāo)和位移待求點坐標(biāo)有關(guān)的函數(shù)，稱為形函數(shù)。在這里，形函數(shù)采用一次多項式的形式，表達(dá)形式為：

式（3）中，D 代表三角形單元的面積，x、y 表示單元任意一點的坐標(biāo)，其他系數(shù)的表示形式為：

式（4）中，“（i，j，k）”的意思是依次將下標(biāo)做輪換：i→j，j→k，k→i。這樣，所有的待定參數(shù)都可以通過節(jié)點位置來表示。

2.2.2 單元應(yīng)變

在確定了單元內(nèi)任意一點的位移之后，通過幾何方程即可確定單元內(nèi)的應(yīng)變：

式（5）中，∈為單元應(yīng)變矩陣，B 為形函數(shù)導(dǎo)數(shù)矩。在三節(jié)點三角形單元中，B 的表示形式為：

2.2.3 單元應(yīng)力

求出單元應(yīng)變以后，利用物理方程可以求出單元內(nèi)的應(yīng)力分布：

在二維平面情況下，結(jié)構(gòu)可以分為平面應(yīng)力問題和平面應(yīng)變問題。對于平面應(yīng)力問題，D 可以表示為：

對于平面應(yīng)變問題，D 可以表示為：

式（8）和式（9）中，E 為彈性模量，μ 為剪切模量。

2.2.4 求解單元剛度矩陣

在求出應(yīng)變和應(yīng)力的表達(dá)方式以后，通過虛位移原理可以得到單元剛度矩陣和單元節(jié)點力列陣。虛位移原理可以表示為：

式（10）中，f 和F 分別表示體積力和面積力。在離散化結(jié)構(gòu)中，單元內(nèi)的虛位移原理可以表示為：

經(jīng)過化簡后，可以得到單元內(nèi)部的有限元方程：

2.3 整體分析

在整體分析的過程中，首先將單元剛度矩陣組裝成整體剛度矩陣，然后施加邊界條件和外加載荷，最后計算求解出節(jié)點位移。

2.3.1 求解整體剛度矩陣

在求出單元剛度矩陣和單元節(jié)點力列陣后，將若干個單元的單元剛度矩陣和節(jié)點力列陣進(jìn)行組裝，即：

式（13）就是結(jié)構(gòu)受力分析的有限元方程，將其簡化書寫為：

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

2.3.2 施加邊界條件和外載荷

邊界條件通常是限制結(jié)構(gòu)在某些位置上的位移，表達(dá)形式為：

外載荷的加載形式有位移加載和力加載兩種形式。兩種施加方式的求解過程有所不同，下面分別推導(dǎo)外載荷形式為力載荷和位移載荷下的有限元方程的求解過程。

（1）外載荷形式為力加載的有限元求解過程

當(dāng)外載荷形式為力載荷時，力載荷的加載形式可以表示為：

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

式（18）中，u1和 P2是已知的參數(shù)，u2和 u3為待求的節(jié)點位移值。此時，可以將有限元方程簡化為：

由式（19）即可求得待求的節(jié)點位移u2和u3。

（2）外載荷形式為位移載荷的有限元求解過程

當(dāng)外載荷形式為位移載荷時，位移載荷的加載形式可以表示為：

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

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

由式（22）即可求得待求的節(jié)點位移u3，然后將求出的節(jié)點位移值帶入到式（21）中即可求得力載荷P2。

2.4 離散體的受力分析

在求出離散體的節(jié)點位移以后，首先通過式（5）求解出每個單元的應(yīng)變分布，然后通過式（7）求解出每個單元的應(yīng)力分布。

3 數(shù)值仿真實例

在數(shù)值算例中，我們運用MATLAB 軟件編寫出力載荷和位移載荷下的有限元程序，在MATLAB 軟件中完成了力載荷和位移載荷下的有限元仿真，并且得到相關(guān)的數(shù)值結(jié)果。

3.1 外載荷形式為力加載的數(shù)值實例

如圖4 所示為二維平面的矩形板。矩形板的長為3m，寬為 1m，彈性模量為 E=1×107Pa，泊松比 μ=0.3，在矩形板右上方的頂點處施加的力載荷大小為F=1×105N。使用三節(jié)點三角形單元，并采用平面應(yīng)力假設(shè)，求解矩形板每個節(jié)點上的位移。

圖4 力載荷作用下的矩形板及其有限元幾何模型

受力載荷作用下的有限元分析步驟如下所示：

（1）將矩形板進(jìn)行離散化處理；

（2）使用式（12）求解出每一個單元的單元剛度矩陣；

（3）使用式（13）組裝整體剛度矩陣；

（4）處理邊界條件和力載荷

矩形板的邊界條件為：

u7=v7=u8=v8=0；

外加的節(jié)點力載荷為：

Py1=-1×104N；

（5）求解未知的節(jié)點位移

運用式（18）和式（19），求出節(jié)點位移列陣為：

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）求解單元內(nèi)的應(yīng)力狀態(tài)

將節(jié)點位移帶入到每一個單元內(nèi)，可求出力載荷情況下單元1～6 中的應(yīng)力狀態(tài)，具體數(shù)值如表1 所示。

3.2 外載荷形式為位移加載的數(shù)值實例

如圖5 所示為一個二維平面的方形板。矩形板的長為3m，寬為 1m，彈性模量為 E=1×107Pa，泊松比 μ=0.3，在矩形板右上方的頂點處施加的位移載荷大小為u=-0.0329m。使用三節(jié)點三角形單元，并采用平面應(yīng)力假設(shè)，求解矩形板每個節(jié)點上的位移和每個單元的應(yīng)力值。

圖5 位移載荷作用下的矩形板及其有限元幾何模型

受位移載荷作用下的有限元分析步驟如下所示：

（1）將矩形板進(jìn)行離散化處理；

（2）使用式（12）求解出每一個單元的單元剛度矩陣；

（3）使用式（13）組裝整體剛度矩陣；

表1 力載荷下各單元的應(yīng)力值（Pa）

表2 位移載荷下各單元的應(yīng)力值（Pa）

（4）處理邊界條件和位移載荷

矩形板的邊界條件為：

u7=v7=u8=v8=0；

外加的節(jié)點位移載荷為：

v1=-0.0329m；

（5）求解未知的節(jié)點位移和外加力載荷

使用式（21）和式（22）求解未知的位移，求出的節(jié)點位移列陣為：

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）求解各單元的應(yīng)力情況

將節(jié)點位移帶入到每一個單元內(nèi)，可求出位移載荷情況下單元1～6 中的應(yīng)力狀態(tài)，具體數(shù)值如表2 所示。

3.3 兩種外加載荷的仿真結(jié)果比較

使用力載荷加載計算出的位移結(jié)果和使用位移載荷計算出的位移結(jié)果幾乎一致，兩種外加載荷計算出來的應(yīng)力結(jié)果也大致相同。說明使用兩種不同的外加載荷進(jìn)行有限元分析得到的數(shù)值結(jié)果是一致的，驗證了運用位移載荷進(jìn)行有限元分析的正確性。

4 結(jié)論

本文分析了外加載荷分別為力載荷和位移載荷情況下的有限元原理，并且使用MATLAB 軟件編寫了兩種外加載荷在二維平面應(yīng)力問題下的有限元程序，驗證了運用力載荷和位移載荷對結(jié)構(gòu)進(jìn)行有限元分析時仿真結(jié)果的一致性。為結(jié)構(gòu)的有限元分析擴展了思路。當(dāng)結(jié)構(gòu)的外加力載荷不容易得到時，通過測量位移載荷的大小也可以對結(jié)構(gòu)進(jìn)行有限元分析。