第三章 有限元法

第3章有限元法

FiniteElementMethod

(FEM)第3章有限元法內(nèi)容簡介

有限元法是結(jié)構(gòu)分析的一種數(shù)值計算方法。它在20世紀(jì)50年代初期隨著計算機(jī)的發(fā)展應(yīng)運(yùn)而生。

這一方法的理論基礎(chǔ)牢靠，物理概念清晰，解題效率高，適應(yīng)性強(qiáng)，目前已成為機(jī)械產(chǎn)品動、靜、熱特性分析的重要手段，它的程序包是機(jī)械產(chǎn)品計算機(jī)輔助設(shè)計方法庫中不可缺少的內(nèi)容之一。

本章介紹了如下內(nèi)容:

■

有限元法的概述

■

平面剛架的有限元法

■

彈性力學(xué)平面問題的有限元法3.1概述

在工程分析和科學(xué)研究中，常常會遇到大量的由常微分方程、偏微分方程及相應(yīng)的邊界條件描述的場問題，如位移場、應(yīng)力場和溫度場等問題。目前求解這類場問題的方法主要有兩種：

●

用解析法求得精確解；

●

用數(shù)值解法求其近似解。其中，能用解析法求出精確解的只能是方程性質(zhì)比較簡單且?guī)缀芜吔缦喈?dāng)規(guī)則的少數(shù)問題。

而對于絕大多數(shù)問題，則很少能得出解析解。這就需要研究它的數(shù)值解法，以求出近似解。

目前，工程中實(shí)用的數(shù)值解法主要有三種：

有限差分法

有限元法

邊界元法其中，以有限元法通用性最好，解題效率高，工程應(yīng)用最廣。目前它已成為機(jī)械產(chǎn)品動、靜、熱特性分析的重要手段，它的程序包是機(jī)械產(chǎn)品計算機(jī)輔助設(shè)計方法庫中不可缺少的內(nèi)容之一。“有限元法

”的基本思想早在20世紀(jì)40年代初期就有人提出，但真正用于工程中則是電子計算機(jī)出現(xiàn)以后?！坝邢拊?/p>

”這一名稱是1960年美國的克拉夫（Clough，R.W.）在一篇題為“平面應(yīng)力分析的有限元法”論文中首先使用。此后，有限元法的應(yīng)用得到蓬勃發(fā)展。到20世紀(jì)80年代初期國際上較大型的結(jié)構(gòu)分析有限元通用程序多達(dá)幾百種，從而為工程應(yīng)用提供了方便條件。由于有限元通用程序使用方便，計算精度高，其計算結(jié)果已成為各類工業(yè)產(chǎn)品設(shè)計和性能分析的可靠依據(jù)。

3.1概述3.1.2桿系單元定義桿系結(jié)構(gòu)中的桿件、梁、柱等稱為桿系單元。連接的點(diǎn)稱為節(jié)點(diǎn)。桿系單元為一維單元。結(jié)構(gòu)離散一般原則：桿系的交叉點(diǎn)、邊界點(diǎn)、集中力作用點(diǎn)、桿件截面尺寸突變處等都應(yīng)該設(shè)置節(jié)點(diǎn)，節(jié)點(diǎn)之間的桿件即構(gòu)成單元。F節(jié)點(diǎn)1節(jié)點(diǎn)2單元①節(jié)點(diǎn)3節(jié)點(diǎn)2單元②有限元法的分析過程可概括如下：●

連續(xù)體離散化●

單元分析●

整體分析●

確定約束條件●

有限元方程求解●

結(jié)果分析與討論

1.連續(xù)體離散化

連續(xù)體：是指所求解的對象（如物體或結(jié)構(gòu)）。

離散化（劃分網(wǎng)格或網(wǎng)絡(luò)化）：是將所求解的對象劃分為有限個具有規(guī)則形狀的微小塊體，把每個微小塊體稱為單元，相鄰兩個單元之間只通過若干點(diǎn)互相連接，每個連接點(diǎn)稱為節(jié)點(diǎn)。相鄰單元只在節(jié)點(diǎn)處連接，載荷也只通過節(jié)點(diǎn)在各單元之間傳遞，這些有限個單元的集合體，即原來的連續(xù)體。

單元劃分后，給每個單元及節(jié)點(diǎn)進(jìn)行編號；

選定坐標(biāo)系，計算各個節(jié)點(diǎn)坐標(biāo)；

確定各個單元的形態(tài)和性態(tài)參數(shù)以及邊界條件等。

圖3-1所示為將一懸臂梁建立有限元分析模型的例子。圖中將該懸臂梁劃分為許多三角形單元；

三角形單元的三個頂點(diǎn)都是節(jié)點(diǎn)。圖3-1懸臂梁及其有限元模型

2.單元分析

連續(xù)體離散化后，即可對單元體進(jìn)行特性分析，簡稱為單元分析。

單元分析工作主要有兩項：

(1)選擇單元位移模式(位移函數(shù))

用節(jié)點(diǎn)位移來表示單元體內(nèi)任一點(diǎn)的位移、應(yīng)變和應(yīng)力，就需搞清各單元中的位移分布。

一般是假定單元位移是坐標(biāo)的某種簡單函數(shù)，用其模擬內(nèi)位移的分布規(guī)律，這種函數(shù)就稱為位移模式或位移函數(shù)。通常采用的函數(shù)形式多為多項式。

根據(jù)所選定的位移模式，就可以導(dǎo)出用節(jié)點(diǎn)位移來表示單元體內(nèi)任一點(diǎn)位移的關(guān)系式。

進(jìn)行單元力學(xué)特性分析，將作用在單元上的所有力（表面力、體積力、集中力）等效地移置為節(jié)點(diǎn)載荷；

采用有關(guān)的力學(xué)原理建立單元的平衡方程，求得單元內(nèi)節(jié)點(diǎn)位移與節(jié)點(diǎn)力之間的關(guān)系矩陣

單元剛度矩陣。

(2)

分析單元的特性，建立單元剛度矩陣

3.

整體分析

把各個單元的剛度矩陣集成為總體剛度矩陣，以及將各單元的節(jié)點(diǎn)力向量集成總的力向量，求得整體平衡方程。

集成過程所依據(jù)的原理是節(jié)點(diǎn)變形協(xié)調(diào)條件和平衡條件。

4.

確定約束條件

由上述所形成的整體平衡方程是一組線性代數(shù)方程，在求解之前，必修根據(jù)具體情況分析，確定求解對象問題的邊界約束條件，并對這些方程進(jìn)行適當(dāng)修正。5.

有限元方程求解通過求解整體平衡方程，即可求得各節(jié)點(diǎn)的位移，進(jìn)而根據(jù)位移可計算單元的應(yīng)力及應(yīng)變。6.

結(jié)果分析與討論F節(jié)點(diǎn)2節(jié)點(diǎn)1單元②節(jié)點(diǎn)1節(jié)點(diǎn)2單元①節(jié)點(diǎn)2節(jié)點(diǎn)1節(jié)點(diǎn)1節(jié)點(diǎn)2節(jié)點(diǎn)1

網(wǎng)架桿件節(jié)點(diǎn)位移單元剛度矩陣總剛度矩陣總剛度方程節(jié)點(diǎn)位移值桿件內(nèi)力單元內(nèi)力與節(jié)點(diǎn)位移間關(guān)系引入邊界條件節(jié)點(diǎn)平衡及變形協(xié)調(diào)條件基本單元基本未知量總疊加cA疊加3.2平面剛架的有限元法桿系單元分類桁架單元：桁架中的桿件剛架單元：剛架中的桿件區(qū)別：桁架節(jié)點(diǎn)：鉸節(jié)點(diǎn)傳遞力！剛架節(jié)點(diǎn)：剛節(jié)點(diǎn)

傳遞力和力矩！

桿系單元的有限元分析與平面問題和空間問題比較，基本流程完全相同；具體計算細(xì)節(jié)需要按照桿系單元的特性來進(jìn)行。剛架的有限元分析平面剛架兩個坐標(biāo)系：局部坐標(biāo)系整體坐標(biāo)系

梁在橫向外載荷（可以是集中力或分布力或力矩等）作用下產(chǎn)生彎曲變形，在水平載荷作用下產(chǎn)生線位移。

對于該平面簡支梁問題：

梁上任一點(diǎn)受有三個力的作用：水平力Fx，

剪切力Fy，

和彎矩Mz。

相應(yīng)的位移為：水平線位移u，

撓度v，

和轉(zhuǎn)角z

。由上圖可見：

水平線位移和水平力向右為正，

撓度和剪切力向上為正，

轉(zhuǎn)角和彎矩逆時針方向為正。

通常規(guī)定：

水平線位移和水平力向右為正，

撓度和剪切力向上為正，

轉(zhuǎn)角和彎矩逆時針方向為正。

為使問題簡化，可把圖示的梁看作是一個梁單元。如圖所示，當(dāng)令左支承點(diǎn)為節(jié)點(diǎn)1

，右支承點(diǎn)為節(jié)點(diǎn)2

時，則該單元的節(jié)點(diǎn)位移和節(jié)點(diǎn)力可以分別表示為：稱為單元的節(jié)點(diǎn)位移列陣。稱為單元的節(jié)點(diǎn)力列陣；若{f}為外載荷，則稱為載荷列陣。

（5-1）（5-2）寫成矩陣形式為顯然，梁的節(jié)點(diǎn)力和節(jié)點(diǎn)位移是有聯(lián)系的。在彈性小變形范圍內(nèi)，這種關(guān)系是線性的，可用下式表示

或

上式稱為單元有限元方程，或稱為單元剛度方程，它代表了單元的載荷與位移之間（或力與變形之間）的聯(lián)系；式中，[K](e)稱為單元剛度矩陣，它是單元的特性矩陣。

對于所示的平面梁單元問題，利用材料力學(xué)中的桿件受力與變形間的關(guān)系及疊加原理，可以直接計算出單元剛度矩陣[K](e)中的各系數(shù)kst(s,t=i,j)的數(shù)值，具體方法如下：

(1)假設(shè)u1=１，其余位移分量均為零，即v1=1=u2=v2=2=0，此時梁單元如圖(a)所示，由梁的變形公式得圖(a)伸縮：撓度：轉(zhuǎn)角：由上述三式可以解得根據(jù)靜力平衡條件

解得

(2)同理，設(shè)v1=１，其余位移分量均為零，即u1=1=u2=v2=2=0,此時梁單元如圖(b)所示，由梁的變形公式得伸縮：撓度：轉(zhuǎn)角：圖(b)由上述三式可以解得利用靜力平衡條件

解得

(3)同理，設(shè)，其余位移分量均為零，即此時梁單元如圖(c)所示，由梁的變形公式得伸縮：撓度：轉(zhuǎn)角：由上述三式可以解得圖(c)利用靜力平衡條件，可得

解得

剩余三種情況，仿此可推出。最后可以得到平面彎曲梁元的單元剛度矩陣為這樣，便求得單元剛度矩陣中前三列剛度系數(shù)?？梢钥闯?，[k](e)為對稱矩陣。XYαyxuvUV整體坐標(biāo)系OXY局部坐標(biāo)系Oxy坐標(biāo)轉(zhuǎn)換從整體坐標(biāo)到局部坐標(biāo)的坐標(biāo)變換矩陣[T]單元剛度矩陣具有下列特點(diǎn)：矩陣具有對稱性矩陣具有奇異性矩陣具有分塊性總剛矩陣具有下列特點(diǎn)：矩陣具有對稱性，計算時不必將所有元素列出，只列出上三角或下三角即可。矩陣具有奇異性矩陣具有稀疏性。XYO1234(1)(2)(3)(4)總剛矩陣中邊界條件的處理方法

未引入邊界條件前，總剛矩陣[K]是奇異的，不能進(jìn)行求解。引入結(jié)構(gòu)邊界條件消除剛體位移后，總剛矩陣為正定矩陣。

位移為零

彈性約束

指定位移

處理方法約束條件的處理位移分量的值為零：置1賦0法

位移分量等于一個已知的非零值：置大數(shù)法非節(jié)點(diǎn)載荷的處理（1）載荷移置原理（2）固定端反力和反力矩的計算平面桁架的有限元分析在整體坐標(biāo)系下的單元剛度矩陣為：n個節(jié)點(diǎn)類型局部坐標(biāo)系下的單元剛度矩陣[k](e)整體坐標(biāo)系下的單元剛度矩陣[K](e)坐標(biāo)轉(zhuǎn)換矩陣[T](e)總體剛度矩陣[K]剛架6*66*66*63n*3n桁架2*24*42*42n*2n【例】平面桁架結(jié)構(gòu)中，某單元局部編碼依次對應(yīng)的總體編碼為8，6，則單元剛度矩陣中的元素k34應(yīng)放入總體剛度矩陣[K]中的行列2003200420052006200790度時，-BAC-EDF先看行后看列200820092010結(jié)構(gòu)離散化位移函數(shù)單元剛度矩陣載荷移植結(jié)構(gòu)總體剛度方程位移邊界條件的處理應(yīng)力計算公式推廣3.3平面問題的有限單元法彈性力學(xué)平面問題的有限元位移法

有限元位移法求解彈性力學(xué)平面問題：離散研究對象，對它進(jìn)行編號，然后建立線性代數(shù)方程組[K]{δ}={P}需要建立剛度矩陣[K]，建立載荷列陣{P}，引入約束邊界條件，才能解方程得出位移{δ},求得應(yīng)力{σ}。彈性力學(xué)平面問題的有限元位移法求解過程與平面桿件系統(tǒng)的求解過程相似，實(shí)際上平面桿件系統(tǒng)的求解過程運(yùn)用了材料力學(xué)中桿件和梁的己知變形關(guān)系，但彈性力學(xué)平面問題形狀復(fù)雜，受力情況多變，沒有可利用的己知變形關(guān)系，在建立剛度矩陣時必須對離散結(jié)構(gòu)物所用的單元位移(變形)進(jìn)行假設(shè)，因此用最簡單的線性關(guān)系假設(shè)，離散平面問題的單元，三角形單元。

平面應(yīng)力問題：只在平面內(nèi)有應(yīng)力，與該面垂直方向的應(yīng)力可忽略，例如薄板拉壓問題。平面應(yīng)變問題：只在平面內(nèi)有應(yīng)變，與該面垂直方向的應(yīng)變可以忽略。對于平面應(yīng)力問題，彈性矩陣[D]為對于平面應(yīng)變問題，只要將式中的E換成E/1-2，換成/1-，即得到其彈性矩陣簡述平面應(yīng)力問題和平面應(yīng)變問題的區(qū)別？(1)應(yīng)力狀態(tài)不同：平面應(yīng)力問題中平板的厚度與長度、高度相比尺寸小很多，所受的載荷都在平面內(nèi)并沿厚度方向均勻分布，可以認(rèn)為厚度方向的應(yīng)力為零。平面應(yīng)變問題中由于z向尺寸大，該方向上的變形是被約束住的，沿z向的應(yīng)變?yōu)榱?。?）彈性矩陣不同：將平面應(yīng)力問題中的E,就成為平面應(yīng)變問題的彈性矩陣。在有限元分析中，為什么要采用半帶存儲？(1)單元尺寸越小，單元數(shù)越多，分析計算精度越高單元數(shù)越多，總剛度矩陣階數(shù)越高，所需計算機(jī)的內(nèi)存量和計算量越大。(2)總體剛度矩陣具有對稱性、稀疏性以及非零元素帶形分布規(guī)律。(3)只存儲主對角線元素以及上（或下）三角矩陣中寬為NB的斜帶區(qū)內(nèi)的元素，可以大大減小所需內(nèi)存量。

彈性體和有限元計算模型X132u1v1u2v2v(x,y).u(x,y)u3v3(x,y)YO

彈性體和有限元計算模型(1)(2)(3)(4)123456圖3-6a平面問題的有限單元法返回圖3-6b半帶寬B=（相鄰節(jié)點(diǎn)號的最大差值D+1）*2平面問題的有限單元法200320042005200620072008200920105.2.2

單元特性的推導(dǎo)方法

單元剛度矩陣的推導(dǎo)是有限元分析的基本步驟之一。目前，建立單元剛度矩陣的方法主要有以下四種：直接剛度法虛功原理法能量變分法加權(quán)殘數(shù)法

梁在橫向外載荷（可以是集中力或分布力或力矩等）作用下產(chǎn)生彎曲變形，在水平載荷作用下產(chǎn)生線位移。

對于該平面簡支梁問題：

梁上任一點(diǎn)受有三個力的作用：水平力Fx，

剪切力Fy，

和彎矩Mz。

相應(yīng)的位移為：水平線位移u，

撓度v，

和轉(zhuǎn)角z

。由上圖可見：

水平線位移和水平力向右為正，

撓度和剪切力向上為正，

轉(zhuǎn)角和彎矩逆時針方向為正。

通常規(guī)定：

為使問題簡化，可把圖示的梁看作是一個梁單元。如圖5-4所示，當(dāng)令左支承點(diǎn)為節(jié)點(diǎn)i

，右支承點(diǎn)為節(jié)點(diǎn)

j

時，則該單元的節(jié)點(diǎn)位移和節(jié)點(diǎn)力可以分別表示為：稱為單元的節(jié)點(diǎn)位移列陣。稱為單元的節(jié)點(diǎn)力列陣；若{F}為外載荷，則稱為載荷列陣。

（5-1）（5-2）寫成矩陣形式為顯然，梁的節(jié)點(diǎn)力和節(jié)點(diǎn)位移是有聯(lián)系的。在彈性小變形范圍內(nèi)，這種關(guān)系是線性的，可用下式表示

或（5-3b）（5-3a）

上式(5-3b)稱為單元有限元方程，或稱為單元剛度方程，它代表了單元的載荷與位移之間（或力與變形之間）的聯(lián)系；式中，[K](e)稱為單元剛度矩陣，它是單元的特性矩陣。

對于圖5-4所示的平面梁單元問題，利用材料力學(xué)中的桿件受力與變形間的關(guān)系及疊加原理，可以直接計算出單元剛度矩陣[K](e)中的各系數(shù)kst(s,t=i,j)的數(shù)值，具體方法如下：

(1)假設(shè)ui=１，其余位移分量均為零，即vi=iz=uj=vj=zj=0，此時梁單元如圖5-5(a)所示，由梁的變形公式得圖5-5(a)伸縮：撓度：轉(zhuǎn)角：由上述三式可以解得根據(jù)靜力平衡條件

由式(5-3a)解得

(2)同理，設(shè)vi=１，其余位移分量均為零，即ui=iz=uj=vj=zj=0,此時梁單元如圖5-5(b)所示，由梁的變形公式得伸縮：撓度：轉(zhuǎn)角：圖5-5(b)由上述三式可以解得利用靜力平衡條件由式(5-3a)解得

(3)同理，設(shè)，其余位移分量均為零，即此時梁單元如圖5-5(c)所示，由梁的變形公式得伸縮：撓度：轉(zhuǎn)角：由上述三式可以解得圖5-5(c)利用靜力平衡條件，可得

由式(5-3a)解得

剩余三種情況，仿此可推出。最后可以得到平面彎曲梁元的單元剛度矩陣為這樣，便求得式（5-3a）單元剛度矩陣中前三列剛度系數(shù)?？梢钥闯觯琜K](e)為對稱矩陣。（5-4）(1)設(shè)定位移函數(shù)

按照有限元法的基本思想：首先需設(shè)定一種函數(shù)來近似表達(dá)單元內(nèi)部的實(shí)際位移分布，稱為位移函數(shù)，或位移模式。

三節(jié)點(diǎn)三角形單元有6個自由度，可以確定6個待定系數(shù)，故三角形單元的位移函數(shù)為（5-5）式(5-5)為線性多項式，稱為線性位移函數(shù)，相應(yīng)的單元稱為線性單元。上式(5-5)也可用矩陣形式表示，即

式中，wkr5avd為單元內(nèi)任意點(diǎn)的位移列陣。

（5-6）

由于節(jié)點(diǎn)i、j、m

在單元上，它們的位移自然也就滿足位移函數(shù)式(5-5)。設(shè)三個節(jié)點(diǎn)的位移值分別為(ui,vi)、(uj,vj)、(um,vm)，將節(jié)點(diǎn)位移和節(jié)點(diǎn)坐標(biāo)代入式(5-5)，得（5-7）式中（5-8）

由上可知，共有6個方程，可以求出6個待定系數(shù)。解方程，求得各待定系數(shù)和節(jié)點(diǎn)位移之間的表達(dá)式為

為三角形單元的面積。其中：

（5-9）將式(5-7)及式(5-8)、式(5-9)代入式(5-6)中，得到（5-10）（5-11）式中，矩陣[N]稱為單元的形函數(shù)矩陣；

為單元節(jié)點(diǎn)位移列陣。其中，為單元的形函數(shù)，它們反映單元內(nèi)位移的分布形態(tài)，是x,y坐標(biāo)的連續(xù)函數(shù)，且有（5-12）式(5-10)又可以寫成（5-13）上式清楚地表示了單元內(nèi)任意點(diǎn)位移可由節(jié)點(diǎn)位移插值求出。

(2)

利用幾何方程由位移函數(shù)求應(yīng)變

根據(jù)彈性力學(xué)的幾何方程

，

線應(yīng)變

剪切應(yīng)變

則應(yīng)變列陣可以寫成式中，[B]

稱為單元應(yīng)變矩陣，它是僅與單元幾何尺寸有關(guān)的常量矩陣，即（5-14）（5-15）

上述方程(5-14)稱為單元應(yīng)變方程，它的意義在于：

單元內(nèi)任意點(diǎn)的應(yīng)變分量亦可用基本未知量即節(jié)點(diǎn)位移分量來表示。

則有如下單元應(yīng)力方程由式(5-19)可求單元內(nèi)任意點(diǎn)的應(yīng)力分量，它也可用基本未知量即節(jié)點(diǎn)位移分量來表示。（5-18）（5-19）(4)由虛功原理求單元剛度矩陣

根據(jù)虛功原理，當(dāng)彈性結(jié)構(gòu)受到外載荷作用處于平衡狀態(tài)時，在任意給出的微小的虛位移上，外力在虛位移上所做的虛功AF等于結(jié)構(gòu)內(nèi)應(yīng)力在虛應(yīng)變上所存儲的虛變形勢能A，即設(shè)處于平衡狀態(tài)的彈性結(jié)構(gòu)內(nèi)任一單元發(fā)生一個微小的虛位移，則單元各節(jié)點(diǎn)的虛位移為（5-21）（5-22）（5-20）則單元內(nèi)部必定產(chǎn)生相應(yīng)的虛應(yīng)變，故單元內(nèi)任一點(diǎn)的虛應(yīng)變?yōu)轱@然，虛應(yīng)變和虛位移之間關(guān)系為設(shè)節(jié)點(diǎn)力為

則外力虛功為

（5-25）（5-23）（5-24）單元內(nèi)的虛變形勢能為

根據(jù)