吳 敏

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

1 概述

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

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

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

2 有限元分析步驟

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

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

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

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

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

2.2 單元分析

將結(jié)構(gòu)經(jīng)過(guò)離散化處理后，結(jié)構(gòu)由若干個(gè)離散的單元表示。在單元分析的過(guò)程中，將用單元的節(jié)點(diǎn)位移分別表示單元內(nèi)部位移、單元應(yīng)變和單元應(yīng)力，并且求出單元?jiǎng)偠染仃嚒?/p>

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

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

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

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

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

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

2.2.2 單元應(yīng)變

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

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

2.2.3 單元應(yīng)力

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

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

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

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

2.2.4 求解單元?jiǎng)偠染仃?/p>

在求出應(yīng)變和應(yīng)力的表達(dá)方式以后，通過(guò)虛位移原理可以得到單元?jiǎng)偠染仃嚭蛦卧?jié)點(diǎn)力列陣。虛位移原理可以表示為：

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

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

2.3 整體分析

在整體分析的過(guò)程中，首先將單元?jiǎng)偠染仃嚱M裝成整體剛度矩陣，然后施加邊界條件和外加載荷，最后計(jì)算求解出節(jié)點(diǎn)位移。

2.3.1 求解整體剛度矩陣

在求出單元?jiǎng)偠染仃嚭蛦卧?jié)點(diǎn)力列陣后，將若干個(gè)單元的單元?jiǎng)偠染仃嚭凸?jié)點(diǎn)力列陣進(jìn)行組裝，即：

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

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

2.3.2 施加邊界條件和外載荷

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

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

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

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

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

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

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

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

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

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

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

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

2.4 離散體的受力分析

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

3 數(shù)值仿真實(shí)例

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

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

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

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

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

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

（2）使用式（12）求解出每一個(gè)單元的單元?jiǎng)偠染仃嚕?/p>

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

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

矩形板的邊界條件為：

u7=v7=u8=v8=0；

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

Py1=-1×104N；

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

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

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é)點(diǎn)位移帶入到每一個(gè)單元內(nèi)，可求出力載荷情況下單元1～6 中的應(yīng)力狀態(tài)，具體數(shù)值如表1 所示。

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

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

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

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

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

（2）使用式（12）求解出每一個(gè)單元的單元?jiǎng)偠染仃嚕?/p>

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

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

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

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

矩形板的邊界條件為：

u7=v7=u8=v8=0；

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

v1=-0.0329m；

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

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

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），求得施加位移載荷時(shí)等效的外加載荷，其值為：

Py1=-0.993×103N；

（6）求解各單元的應(yīng)力情況

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

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

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

4 結(jié)論

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