有限元分析基礎(chǔ)

FundamentalsofFiniteElementAnalysis

唐文亭

133892729582011.9

內(nèi)容結(jié)構(gòu)第一章概述第六章軸對稱問題的有限單元法第七章軟件介紹及計(jì)算示例第五章空間問題的有限單元法第四章等參元第三章桿系結(jié)構(gòu)靜力分析的有限單元法第二章彈性力學(xué)平面問題有限元法第二章彈性力學(xué)平面問題有限元法有限元法的基本思想及優(yōu)越性在應(yīng)用有限元法時，我們首先將一個連續(xù)的彈性體看作由許多尺寸有限的小單元---有限元組成。

這就是所謂區(qū)域劃分，在數(shù)學(xué)上稱為“離散化”。2.根據(jù)計(jì)算對象的簡化模型，單元的形狀取成平面三角形或四邊形，四面體或六面體等。單元與單元之間，通過若干個稱為“節(jié)點(diǎn)”的點(diǎn)鉸接相連，由此組合成整體。3.以一個個小單元為計(jì)算單位，首先進(jìn)行單元分析，然后把它們組裝起來，進(jìn)行整體分析，最后求出結(jié)構(gòu)的近似解。

這種把復(fù)雜結(jié)構(gòu)看成有限個單元組成的整體，就是有限元法的基本思想。

有限元法是從基于能量的變分法發(fā)展而來的。如應(yīng)用最小勢能原理的雷利---里茲法，當(dāng)按位移求解時，它首先要尋找一個滿足整個彈性體幾何邊界條件的位移函數(shù)，這對工程實(shí)際問題往往有困難。

而用有限元法時，將結(jié)構(gòu)進(jìn)行離散，從一個個單元入手，只要假設(shè)單元上的分片插值函數(shù)，然后綜合起來，代替整個域上的位移函數(shù)，這就使問題大為簡便和靈活。

因此，有限元法是以變分原理和分片插值為基礎(chǔ)的。得到的是近似解。（1）有限元法直接在力學(xué)模型上進(jìn)行離散化（剖分），

物理概念清晰，明白易懂。（2）有限元法有較好的適應(yīng)性。對于簡單問題和復(fù)雜問題基本上同等處理。（3）有限元法的各個計(jì)算步驟，如單元分析，總體分析和方程解算等都較易標(biāo)準(zhǔn)化和程式化,有一套比較固定的分析順序，目前已發(fā)展成各種通用程序，便于掌握和使用。有限元法應(yīng)用于應(yīng)力分析，按所選取的未知量不同可分為三類：

(1)位移法--取節(jié)點(diǎn)位移作為基本未知量；

(2)力法--取節(jié)點(diǎn)力作為基本未知量；

(3)混合法--取一部分節(jié)點(diǎn)位移和一部分節(jié)點(diǎn)力

作為基本未知量。

在推導(dǎo)有限元方程時，主要有兩種方法：

直接法（如直接剛度法）；

變分法（如固體力學(xué)中的最小勢能原理和最小余能原理）把問題歸結(jié)為求泛函的極值問題。

作為初步介紹，我們將以直接剛度法來討論彈性力學(xué)平面問題中的有限元法概念。有限元模型是真實(shí)系統(tǒng)理想化的數(shù)學(xué)抽象。真實(shí)系統(tǒng)有限元模型節(jié)點(diǎn)和單元節(jié)點(diǎn)：空間中的坐標(biāo)位置，具有一定自由度和存在相互物理作用單元:一組節(jié)點(diǎn)自由度間相互作用的數(shù)值矩陣描述(稱為剛度或系數(shù)矩陣)單元有線、面或?qū)嶓w以及二維或三維的單元等種類。載荷10

有限元模型由一些簡單形狀的單元組成，單元之間通過節(jié)點(diǎn)連接，并承受一定載荷。1.每個單元的特性是通過一些線性方程式來描述的。2.作為一個整體，單元形成了整體結(jié)構(gòu)的數(shù)學(xué)模型。3.盡管梯子的有限元模型低于100個方程（即“自由度”），然而在今天一個小的ANSYS分析就可能有5000個未知量，矩陣可能有25,000，000個剛度系數(shù)。節(jié)點(diǎn)自由度是隨連接該節(jié)點(diǎn)單元類型變化的JIIJJKLILKIPOMNKJIL三維桿單元(鉸接)UX,UY,UZ三維梁單元二維或軸對稱實(shí)體單元UX,UY三維四邊形殼單元UX,UY,UZ,三維實(shí)體熱單元TEMPJPOMNKJIL三維實(shí)體結(jié)構(gòu)單元ROTX,ROTY,ROTZROTX,ROTY,ROTZUX,UY,UZ,UX,UY,UZ12用有限元法對彈性力學(xué)平面問題進(jìn)行應(yīng)力分析，不僅具有實(shí)際意義，而且?guī)в幸欢ǖ牡湫托浴Ｍㄟ^它可以看到:

(1)一般情況下處理問題的方法

(2)有限元法的特點(diǎn)

(3)使用中的應(yīng)注意的問題

可為今后進(jìn)一步的深入研究打下基礎(chǔ)。

當(dāng)取節(jié)點(diǎn)位移為基本未知量時，有限元法的解題步驟歸納如下：

下面我們就按上述順序介紹。應(yīng)力計(jì)算方程求解整體分析區(qū)域剖分

單元分析

嚴(yán)格地說，任何彈性體都是處于三維受力狀態(tài)，因而都是空間問題，但是在一定條件下，許多空間問題都可以簡化成平面問題。平面問題可以分為兩類：平面應(yīng)力問題和平面應(yīng)變問題。平面問題應(yīng)力狀態(tài)2.1平面應(yīng)力問題和平面應(yīng)變問題平面應(yīng)力問題

如圖所示的深梁結(jié)構(gòu)，其厚度方向的尺寸遠(yuǎn)比其它兩個方向的尺寸小得多，可視為一薄板。它只承受作用在其平面內(nèi)的載荷，且沿厚度方向不變，計(jì)算時以中性面為研究對象。其力學(xué)特點(diǎn)是：平面應(yīng)力問題的應(yīng)力應(yīng)變轉(zhuǎn)換矩陣即彈性矩陣為：。平面應(yīng)力問題平面應(yīng)變問題平面應(yīng)變問題

圖示為一圓形涵洞的橫截面。其長度方向上的尺寸遠(yuǎn)比其它兩個方向上的尺寸大得多，同樣，載荷作用在xy坐標(biāo)面內(nèi)，且沿z軸方向均勻分布。其力學(xué)特點(diǎn)是：但一般情況下:平面應(yīng)變問題的彈性矩陣只需將式(4-1)中的E換成換成

，即可。。無論是平面應(yīng)力問題還是平面應(yīng)變問題的應(yīng)力

與應(yīng)變

之間的關(guān)系均為：

其中：為初應(yīng)變。式中(a)三結(jié)點(diǎn)三角形單元(b)四結(jié)點(diǎn)正方形單元(c)四結(jié)點(diǎn)矩形單元(d)四結(jié)點(diǎn)四邊形單元平面問題單元的主要類型2.2彈性體的剖分

作為用有限元法解決彈性力學(xué)問題的第一步，必須先對彈性體區(qū)域進(jìn)行剖分。對于平面問題來說，最簡單的方法是用直線將彈性體區(qū)域剖分為有限個三角形或四邊形單元。本章將只討論三個節(jié)點(diǎn)的三角形單元。

圖2-1剖分要一直進(jìn)行到彈性區(qū)域的邊界上當(dāng)邊界是直線段時，就取其為三角形單元的一條邊；當(dāng)邊界是曲線時，則在每小段上用相應(yīng)的直線近似地代替曲線而作為三角形單元的一邊，如圖2-1

單元分得越小，結(jié)構(gòu)計(jì)算越精確。因此，應(yīng)當(dāng)在計(jì)算機(jī)容量的允許的范圍內(nèi)，盡可能地提高工程上的精確要求，適當(dāng)?shù)卮_定單元的大小和數(shù)目。

單元的大小和數(shù)目要根據(jù)精度的要求和計(jì)算機(jī)容量來確定。

1）任意一個三角形單元的頂點(diǎn)，必須同時也是其相鄰的三角形的單元的頂點(diǎn)（如圖2-2a)，而不能是其相鄰三角形的內(nèi)點(diǎn)（如圖2-2b)。圖2-2a圖2-2b具體進(jìn)行剖分時，一般應(yīng)注意以下幾點(diǎn)：2）盡可能使同一個三角形單元各邊的長度相差不太大。

此外在三角形單元中最好不要出現(xiàn)鈍角。

因此，在圖2-3a、b兩種剖分式中，雖然都涉及

到同樣的四個頂點(diǎn)，但我們通常都采用a，而不采用b。圖2-3a圖2-3b3）在事先估計(jì)應(yīng)力較為集中、應(yīng)力變化較大的地方，例如孔洞附近以及形狀突變的角點(diǎn)等處，單元應(yīng)分得小一些；在應(yīng)力變化比較平緩的地方，如離

開孔洞一定的距離處，單元可以分得比較大一點(diǎn)，

如圖2-4。圖2-4有時應(yīng)力情況事先無法估計(jì)，可先采用比較均勻的剖分法進(jìn)行一次初算，然后考慮初算的結(jié)果，重新合理剖分，再進(jìn)行第二次計(jì)算，或用光彈性的方法事先對應(yīng)力場作一個大概的了解，再在此基礎(chǔ)上作合理的剖分和計(jì)算，這也是一種常用的方法。4）在厚度或材料常數(shù)有突變的地方，除了應(yīng)把這些部位的單元分得較小，較密一些以外，還必須把突變線作為單元的分界線。也就是說，在一個單元內(nèi)部，只能包含一個厚度和一種材料常數(shù)。

5）當(dāng)整個彈性體區(qū)域在幾何上具有對稱軸，而載荷又對稱于該軸或反對稱于該軸時，則其位移和應(yīng)力也必然具有這種對稱性質(zhì)。為了減少計(jì)算量，只需取其一部分作為求解區(qū)域進(jìn)行單元剖分和計(jì)算即可。

(a)均勻受力板力學(xué)模型(b)力學(xué)模型離散化平面問題有限單元法的計(jì)算力學(xué)模型

作了這樣的剖分之后，再以三角形單元的頂點(diǎn)作為節(jié)點(diǎn)（注意，如果邊界上有集中力，則一般將其作用點(diǎn)選定為節(jié)點(diǎn)），然后對單元和節(jié)點(diǎn)分別進(jìn)行編號。

編號的順序不影響計(jì)算結(jié)果，原則上是可以任意的。但用直接法求解有限元的基本方程時，從壓縮計(jì)算機(jī)存儲量的角度來看，在對節(jié)點(diǎn)編號時應(yīng)注意：單元的兩個相鄰節(jié)點(diǎn)編號之差應(yīng)盡可能地小。因?yàn)檫@個差值就反映在方程組的系數(shù)矩陣（總剛度矩陣）的帶寬上，它直接決定了系數(shù)矩陣元素的存儲數(shù)量。有關(guān)問題，以后還要作詳細(xì)的說明。2.3單元分析

在進(jìn)行了彈性體的剖分后，可任取一單元作為研究對象。設(shè)某三角形單元e的節(jié)點(diǎn)編號為i,j,m,（為了在以后的計(jì)算中使三角形的面積不致為負(fù)值，規(guī)定i,j,m的次序?yàn)槟鏁r針方向。）并設(shè)三個節(jié)點(diǎn)i,j,m在右手坐標(biāo)系的坐標(biāo)值分別為:（xi,yi），（xj，yj），（xm，ym），如圖2-5所示。圖2-5

對于平面問題，三個節(jié)點(diǎn)的位移分別為：單元的節(jié)點(diǎn)位移列陣為：所謂單元分析，就是建立節(jié)點(diǎn)位移{}（基本未知量）和單元內(nèi)任意一點(diǎn)的：位移{f}，單元應(yīng)變{ε}，單元應(yīng)力{σ}單元節(jié)點(diǎn)力{F}e

之間的關(guān)系，使{f}，{ε}，{σ}，{F}e等都用節(jié)點(diǎn)位移{}e來表示。如此，則基本未知量{}e一經(jīng)求得，其它各量皆可隨之而定。

1）節(jié)點(diǎn)位移{}e

和單元內(nèi)任意一點(diǎn)位移{f}關(guān)系

首先，我們要確定三角形單元內(nèi)各點(diǎn)的位移變化規(guī)律。即當(dāng)節(jié)點(diǎn)位移確定時，單元內(nèi)各點(diǎn)的位移應(yīng)如何插值？

設(shè)單元內(nèi)任一點(diǎn)的位移是該點(diǎn)坐標(biāo)（x,y)的線性函數(shù)。對于采用三角形單元的平面問題來說，當(dāng)單元取得足夠小時，取線性位移插值函數(shù)是合理的。即：

式中

是待定常數(shù)，它們可以由單元的邊界條件，即節(jié)點(diǎn)的位移值來確定。

為此，只要將節(jié)點(diǎn)的坐標(biāo)值代入式（a)，就得到節(jié)點(diǎn)的位移值：

用克萊姆法則求解線性方程組（b）,（c），得:

式中：

而:

是三角形單元的面積將式(d)、(e)代入式(a)，即得單元的位移插值函數(shù)：

進(jìn)行整理后得：

若令：

為單元的形函數(shù)。由式(f)可知，單元內(nèi)任意點(diǎn)的位移與單元的節(jié)點(diǎn)位移是通過形函數(shù)來聯(lián)系的，而形函數(shù)則是點(diǎn)的坐標(biāo)的線性函數(shù)。引入式(2-5)后，式(f)可以表示為：寫成矩陣形式就是：

式中

為二階單位矩陣

若

則單元上的位移插值函數(shù)可表示為：

{f}=[N]e{

}e

應(yīng)當(dāng)指出，任意兩個相鄰的三角形單元，如圖中的

i,j,m及p,j,i,它們在i和j點(diǎn)具有相同的位移。

我們已假定了位移分量在每一單元中是坐標(biāo)的線性函數(shù)，則在公共邊ij上，位移也必然是按同樣的線性變化的。

因此，在上述兩個單元中，公共邊ij上各點(diǎn)也都具有相同的位移。

這就保證了相鄰單元在公共邊界上位移的連續(xù)性，也即彈性體在受力變形后，各單元的邊界線上的材料不致產(chǎn)生空隙或重疊的現(xiàn)象。

2)節(jié)點(diǎn)位移{}e和單元應(yīng)變分量{ε}的關(guān)系

從上一節(jié)可知，在取定單元的位移插值函數(shù)以后，只要求得各個節(jié)點(diǎn)的位移值，則每個單元內(nèi)各點(diǎn)的位移（因而也是整個彈性體內(nèi)各點(diǎn)的位移）即可確定。

這一節(jié)，我們將聯(lián)系到平面問題的幾何方程和物理方程，來進(jìn)一步確定單元的應(yīng)變和應(yīng)力。

由(2-4)式，即可得到用節(jié)點(diǎn)位移表示單元任一點(diǎn)的應(yīng)變表達(dá)式：

寫成矩陣形式：簡寫為：

{ε}=[B]{}e

其中：稱為單元的應(yīng)變矩陣。

由于單元的面積△以及各幾何參數(shù)bi,ci,……cm的值，都可以由節(jié)點(diǎn)坐標(biāo)直接確定，而且均為常數(shù)。

因此在每一個單元中，應(yīng)變分量

εx,εy,εxy

都是常量。

故線性位移插值函數(shù)的單元又稱為常應(yīng)變單元。

3）節(jié)點(diǎn)位移{}e和單元應(yīng)力分量{σ}的關(guān)系

由彈性力學(xué)已知，平面應(yīng)力情況下，應(yīng)力與應(yīng)變之間的關(guān)系可表達(dá)為：寫成矩陣形式：或縮寫成：

{σ}=[D]{ε}(2-12)

式中：

[D]—彈性矩陣在平面應(yīng)變情況下，彈性矩陣為：將（2-10）式代入（2-12）式，即可得到用節(jié)點(diǎn)位移表示的單元內(nèi)任意一點(diǎn)應(yīng)力的表達(dá)式：

{σ}=[D][B]{}e

或?qū)懗?

{σ}=[S]{}e

式中:

[S]=[D][B]稱為應(yīng)力矩陣。4）節(jié)點(diǎn)位移{}e和單元節(jié)點(diǎn)力{F}e的關(guān)系---單元剛度矩陣[k]

所謂有限元法的“位移法”，就是把節(jié)點(diǎn)位移作為未知數(shù)。而把單元內(nèi)各點(diǎn)的位移、應(yīng)力和應(yīng)變都表達(dá)成節(jié)點(diǎn)位移的函數(shù)。

這樣，就把全部問題都?xì)w納為求取節(jié)點(diǎn)位移的問題了。

在彈性力學(xué)中，以應(yīng)力形式表示的物體內(nèi)部各部分之間的相互作用力，在有限元法中就相應(yīng)地以節(jié)點(diǎn)力的形式來代替。

所謂節(jié)點(diǎn)力就是單元周圍部分對我們所考慮的那個單元的作用。

設(shè)三角形單元i,j,m三個節(jié)點(diǎn)上的節(jié)點(diǎn)力分量分別為：則單元的節(jié)點(diǎn)力列陣為：與應(yīng)力---應(yīng)變成線性關(guān)系相似，就彈性系統(tǒng)而言，節(jié)點(diǎn)位移和節(jié)點(diǎn)力之間，也同樣保持線性關(guān)系：用矩陣形式表示：或簡寫成：

{F}e={k}e

{}e

(2-16b)

這就是單元節(jié)點(diǎn)位移與節(jié)點(diǎn)力的表達(dá)式。式中：

稱為單元剛度矩陣。關(guān)于它的性質(zhì)，可以討論如下：

①單元剛度矩陣每一列的意義：

（為了討論方便，設(shè)i=1,j=2,m=3）令u1=1,而v1=u2=v2=u3=v3=0，由（2-16a）式可得，由此可見,單元剛度矩陣的第一列元素表示:

當(dāng)節(jié)點(diǎn)1在x方向有單位位移（u1=1）

而其它位移為零（v1=u2=v2=u3=v3=0）時，各節(jié)點(diǎn)上產(chǎn)生的節(jié)點(diǎn)力。

如k61，就表示當(dāng)節(jié)點(diǎn)1在x方向有單位位移時，在節(jié)點(diǎn)3,y方向產(chǎn)生的節(jié)點(diǎn)力。

因單元在這些力作用下處于平衡，所以x方向和y方向的節(jié)點(diǎn)力之和分別為零，其總和亦為零，從而有：

k11+k21+k31+k41+k51+k61=0

對于其它各列也有類似的性質(zhì)，元素kij下標(biāo)j表示產(chǎn)生單位位移的節(jié)點(diǎn)序號和方向，i表示產(chǎn)生節(jié)點(diǎn)力的節(jié)點(diǎn)序號和方向。②單元剛度矩陣主對角線上的元素為正

例如k11,表示節(jié)點(diǎn)1在x方向產(chǎn)生單位位移而其它位移均為零時，必須在節(jié)點(diǎn)1,x方向施加的力。

該力顯然應(yīng)和位移方向一致。因而k11應(yīng)為正值。

③單元剛度矩陣是對稱矩陣它的對稱性是由彈性結(jié)構(gòu)的反力互等定理（第j個單位位移分量引起的第i個節(jié)點(diǎn)力等于第i個單位位移分量引起的第j個節(jié)點(diǎn)力）得到的：

即：

通過單元剛度矩陣[k]e,建立了單元節(jié)點(diǎn)力列陣{F}e與單元節(jié)點(diǎn)位移列陣{}e之間的關(guān)系式（2-16)。

在有限元法中，必須逐一求出每個單元的單元剛度矩陣。這是計(jì)算中必不可少的重要一環(huán)。為了求取單元剛度矩陣，我們應(yīng)用虛功原理推演它的具體表達(dá)式。

任取某一單元，虛功原理可簡述如下：

當(dāng)處于平衡狀態(tài)的單元體發(fā)生約束條件所允許的微小虛位移時，則節(jié)點(diǎn)力在虛位移上所作的虛功等于整個單元體的虛變形能。

(a)實(shí)際力系(b)虛設(shè)位移彈性體虛功原理的應(yīng)用如以{

*}表示節(jié)點(diǎn)虛位移，{ε*}表示與虛位移相應(yīng)的虛應(yīng)變，

即:虛應(yīng)變{ε*}與節(jié)點(diǎn)虛位移{

*}之間的關(guān)系仍然如公式（2-10）所示，即:

{ε*}=[B]{

*}

節(jié)點(diǎn)力{F}e在節(jié)點(diǎn)虛位移{

*}上所作的功為：

整個單元體的虛應(yīng)變能:

(其中t為薄板厚度）按虛功原理，（b）和（c）相等，即:

將（a）式和（2-14）代入上式的右端:

根據(jù)矩陣相乘的逆序法則：

（[B]{

*})T={

*}T[B]T

便有：

由于單元節(jié)點(diǎn)虛位移{

*}是任意給的微小量，故{

*}T可以提到積分符號外面去。于是得：

與(2-16b)式:

比較可得：

由于彈性矩陣[D]完全決定于彈性常數(shù)E和μ，而在選取線性位移插值函數(shù)的條件下，應(yīng)變矩陣[B]中的元素又都是常量，因此可把（h）式中的常量都提到積分號外面來，并注意到:

是三角形ijm的面積，于是:

[k]e=[B]T[D][B]t

△

用矩陣形式表示為：算例：設(shè)某三角形單元ijm如下圖(2-8)所示。

已知：彈性常數(shù)E=1.5,v=0.25;厚度t=0.4

試計(jì)算:

應(yīng)變矩陣[B]

應(yīng)力矩陣[S]

單元剛度矩陣[k]e

1.常數(shù)計(jì)算和面積計(jì)算

面積△=1/2（bicj–bjci)=0.5

2.應(yīng)變矩陣[B]和應(yīng)力矩陣[S]:

3.單元剛度矩陣

[k]e

=[B]T[D][B]△t=[B]T[S]△t

根據(jù)上述條件，試計(jì)算圖(2-9)所示三角形ijm的單元剛度矩陣：

答：由公式(2-2):

可以看出，bi,bj,bm和ci，cj，cm等只與單元的方位有關(guān)，不隨坐標(biāo)的平移而變化。可以利用這一性質(zhì)適當(dāng)選擇原點(diǎn)的位置以簡化單元剛度矩陣的計(jì)算。

即對于右圖。下述兩種坐標(biāo)的選取，所得[k]e一樣。當(dāng)把單元剛度矩陣寫成2x2的分塊子矩陣時，

[k]e可寫成下列形式：式中:如為平面應(yīng)變問題，將(2-19)式中的:

E換成E/(1-μ2)

μ換成μ/(1-μ)

于是得到：

按(2-20a）式或（2-20b）式先依次計(jì)算出分塊矩陣中的四個元素可同樣得到單元剛度矩陣中的36個元素。

在實(shí)際計(jì)算中，每個單元的節(jié)點(diǎn)都有一個編號，稱為總體編碼（或?qū)嶋H編碼）。如圖（2-10）。

如果按實(shí)際編碼進(jìn)行運(yùn)算，就會給編制程序帶來困難；因?yàn)槊總€單元的實(shí)際編碼并無一定規(guī)律。但是如賦予每個單元以局部編碼（1）、（2）、（3）這樣就可使用同一個過程體（子程序）。使之能十分方便地對局部編碼的（1）、（2）、（3）三角形單元進(jìn)行運(yùn)算。①②③④①②③④

總體編碼

123

如單元①

局部編碼

①②③

總體編碼

245

單元②

局部編碼

①②③

其余類推5）等效節(jié)點(diǎn)載荷

作用在彈性體上的載荷不外乎是：

*體力（自重、慣性力）

*面力

*集中力

三種。

用有限元法解題時，既然全部問題都?xì)w結(jié)到節(jié)點(diǎn)來處理，那么，當(dāng)單元上作用有外載荷時，也應(yīng)把它們移置到節(jié)點(diǎn)上來，成為節(jié)點(diǎn)載荷。

這種移置必須按照靜力等效原則進(jìn)行

所謂“靜力等效”，系指原載荷與移置后的節(jié)點(diǎn)載荷，在彈性體的任何虛位移過程中的虛功相等。

當(dāng)插值函數(shù)已經(jīng)確定時，這種移置的結(jié)果是唯一的。

在取線性位移插值時，符合剛體靜力等效原則，即：

載荷與節(jié)點(diǎn)載荷在任一軸上投影之和相等，對任一軸的力矩之和也相等，也就是說，原載荷與節(jié)點(diǎn)載荷將具有相同的主矢和主矩。

通常，我們總將集中力的作用點(diǎn)取為節(jié)點(diǎn)，不需要移置。因此，下面只討論：

體力和面力的移置(1)體力的移置以重力為例來說明這個問題。設(shè)有勻質(zhì)、等厚度、編號為e的三角形單元，三個節(jié)點(diǎn)為i,j,m，重力作用在形心c上，如圖2-11所示。由初等幾何知，

首先，我們求移置到i節(jié)點(diǎn)上的垂直節(jié)點(diǎn)載荷。為了便于計(jì)算虛功，假想該單元在節(jié)點(diǎn)i處沿y方向產(chǎn)生一個單位虛位移，而其它兩點(diǎn)不動；這相當(dāng)于在j點(diǎn)及m點(diǎn)安置了鉸支座，在i點(diǎn)安置了水平連桿支座，如圖2-11。

由于我們在三角形單元中采用的位移插值函數(shù)是線性的，所以任一條直線上各點(diǎn)位移都呈線性變化。

現(xiàn)在m點(diǎn)及j點(diǎn)的位移都等于0，所以在邊上各點(diǎn)位移都等于0；

線上各點(diǎn)的垂直位移也按線性變化，在b點(diǎn)等于0，在i點(diǎn)為1。因此c點(diǎn)的垂直位移將為1/3。

按靜力等效原則，體力載荷W的虛功應(yīng)等于的虛功，即有：

或Yei=-W/3

負(fù)號表示Yei的方向與圖上所畫的方向相反。

用同樣的方法可以得到:

Yej=-W/3,Yem=-W/3

下面來求移置到節(jié)點(diǎn)i上的水平節(jié)點(diǎn)載荷Xei。與前面一樣，在形心c處有W作用，假設(shè)節(jié)點(diǎn)i只沿x方向產(chǎn)生一個單位的虛位移，而其它兩點(diǎn)不動。由圖2-12可知，c點(diǎn)的垂直位移等于0，水平位移等于1/3。按靜力等效原則，有：

故:用同樣的方法可以得到：

Xej=0,Xem=0寫成分量形式:

由此可以得出如下結(jié)論：

對于勻質(zhì)、等厚度的三角形單元，當(dāng)考慮自重時，只需把1/3的重量移置到每個節(jié)點(diǎn)上，就完成了重力載荷的移置，而不必再去列出虛功相等的條件。這也完全符合對剛體的靜力等效原則。

但必須指出：上述結(jié)果是由于我們采用線性位移插值函數(shù)造成的。如果位移插值函數(shù)是非線性的，例如是坐標(biāo)的二次函數(shù)，那就不滿足1/3的關(guān)系，也就不能按簡單的剛體靜力等效原則來處理，而必須用虛功方程來建立普遍的表達(dá)式。(2)面力的移置1）設(shè)等厚度的三角形單元e的三個節(jié)點(diǎn)為i,j,m其邊界ij上受有垂直均勻分布的面力載荷。載荷集度（單位長度上的力）為q，如圖2-13所示。仍然采用線性位移插值函數(shù)方法。則根據(jù)靜力等效原則，將此均勻分布的面力載荷移置到兩側(cè)節(jié)點(diǎn)i,j上時，等效節(jié)點(diǎn)載荷為：

式中l(wèi)為

的邊長，

Rei,Rej仍與原載荷平行，

故此時單元的節(jié)點(diǎn)載荷列陣為：

或:2)若ij邊上受有三角形分布的面力載荷，

在i點(diǎn)上的載荷集度為

,j點(diǎn)上為0，則其單元的節(jié)點(diǎn)載荷列陣是：

或：

3）若ij邊上受有垂直的梯形分布的面力載荷，在i,j點(diǎn)上載荷集度分別為和，如圖2-15所示，則其等效節(jié)點(diǎn)載荷為：或?qū)懗煞至啃问剑?/p>

式中分別是在x,y

方向上的分量，其方向與x,y軸正向一致為正，反之為負(fù)。2.3整體分析1.基本方程總剛度矩陣[K]的形成(1)節(jié)點(diǎn)力的組合

以上我們分析了一個單元的情況，現(xiàn)在進(jìn)而研究單元的組合。為了說明問題，今選用一個包含9個節(jié)點(diǎn)8個單元的平面問題來分析。如圖2-16所示，除1、3、7、9四個節(jié)點(diǎn)外，其余五個節(jié)點(diǎn)均聯(lián)接著四個單元。

對于聯(lián)結(jié)著n個單元的節(jié)點(diǎn)，一個節(jié)點(diǎn)的位移當(dāng)然涉及到n個單元，與節(jié)點(diǎn)位移相應(yīng)的節(jié)點(diǎn)力將是n個單元的綜合效應(yīng)。

如節(jié)點(diǎn)5的位移涉及到(2)、(3)、(6)、(7)單元。與節(jié)點(diǎn)5相應(yīng)的節(jié)點(diǎn)力將與上述四個單元有關(guān)。

首先，我們按(2-16a)式逐個建立單元節(jié)點(diǎn)力和節(jié)點(diǎn)位移的關(guān)系。下面，選寫其中(2)、(3)、(6)、(7)四個單元。對于單元(2)，節(jié)點(diǎn)編號2、5、4對于單元(3)，節(jié)點(diǎn)編號2，6，5對于單元(6),節(jié)點(diǎn)編號4，5，8對于單元(7),節(jié)點(diǎn)編號5，6，8

同樣可寫出(1)、(4)、(5)、(8)單元的節(jié)點(diǎn)力和節(jié)點(diǎn)位移的關(guān)系。(學(xué)習(xí)者可自己完成）。

對于第5個節(jié)點(diǎn)，其X方向的節(jié)點(diǎn)力

即:

U5=U5(2)+

U5(3)+

U5(6)+

U5(7)

將上式代入(a)式中:所以式中的上式可簡寫成：

{F}=[K]{

}

式中[K]稱為總剛度矩陣。

它是一個對稱正定陣(其對角線上各元素為正值，且有元素Kij=Kji)。

總剛度矩陣由單元剛度矩陣?yán)奂佣?，每個單元剛度矩陣對總剛度矩陣都有一定貢獻(xiàn)。

[K]=∑[k]e

假定結(jié)構(gòu)離散化后有n個節(jié)點(diǎn)，每個節(jié)點(diǎn)有兩個方程。因此總剛度矩陣為(2nx2n)的矩陣。

可將單元剛度矩陣用補(bǔ)零的方法由6x6擴(kuò)大到(2nx2n)的方陣(圖中虛點(diǎn)上的元素均為0）

如圖所示，則單元剛度矩陣中各元素在總剛度矩陣中的位置即可確定。

例如將第3單元剛度矩陣中的元素填入總剛度矩陣（亦即將該單元剛度矩陣用補(bǔ)零的方法擴(kuò)大成總剛度矩陣）

在計(jì)算程序編制中，我們的做法是先都按局部編碼，使用同一過程體算出各單元的剛度矩陣；然后轉(zhuǎn)換成總體編碼，最后將相同編碼的元素合并成總剛度矩陣中的元素。轉(zhuǎn)換過程示意如下：局部編碼編碼轉(zhuǎn)換總體編碼相同編碼合并形成總剛度矩陣系數(shù)(2)節(jié)點(diǎn)載荷組合

當(dāng)進(jìn)行單元組合時，除了考慮節(jié)點(diǎn)力的組合外，同時還應(yīng)進(jìn)行節(jié)點(diǎn)載荷的組合。

設(shè)結(jié)構(gòu)上的載荷（如體力、面力、集中力等）均已移置到節(jié)點(diǎn)上，則單元節(jié)點(diǎn)載荷列陣為：YmeXmeYieXieXjeYje圖2-17xy對于聯(lián)結(jié)著兩個以上單元的節(jié)點(diǎn)，把相同方向上的載荷迭加起來，顯然有：這就是該節(jié)點(diǎn)載荷在x和y方向的兩個分量。（Σ表示對環(huán)繞節(jié)點(diǎn)i的單元求和）若各節(jié)點(diǎn)上的Xi和Yi（i=1，2，3…n）均已經(jīng)求出，并按節(jié)點(diǎn)編碼的順序排列起來，就得到彈性體總的節(jié)點(diǎn)載荷列陣：(3)平衡方程

在求得了節(jié)點(diǎn)外力矩陣以后，我們就可以寫出位移法中位移分量必須滿足的平衡條件。在有限元法中，也就是節(jié)點(diǎn)位移必須滿足的節(jié)點(diǎn)平衡條件。

根據(jù)公式（2-21）和（2-22），表示所有節(jié)點(diǎn)內(nèi)力與外力平衡的數(shù)學(xué)表達(dá)式為：通常寫成：

等式右端為總節(jié)點(diǎn)載荷列陣，當(dāng)彈性體上的外載荷確定時，它是已知的；總剛度矩陣[K]由單元剛度矩陣[k]e集合而成。因此，式（2-23a）是一個以{}為未知量，以[K]為系數(shù)的線性代數(shù)方程組，這是有限元法的基本方程。其矩陣展開式為：

順便提一下，為了編程方便，對各矩陣元素的下標(biāo)，均按其所在位置標(biāo)定：即為節(jié)點(diǎn)1，x方向的位移u1，

為節(jié)點(diǎn)1，y方向位移v1,......余此類推。

R1為節(jié)點(diǎn)1，x方向的節(jié)點(diǎn)載荷x1，

R2為節(jié)點(diǎn)1，y方向的節(jié)點(diǎn)載荷y1,......余此類推。2.關(guān)于總剛度矩陣的性質(zhì)和常用的存貯方法：

在有限元法中，結(jié)構(gòu)總剛度矩陣的性質(zhì)和常用的存貯方法：

在有限元法中，結(jié)構(gòu)總剛度矩陣的階數(shù)為節(jié)點(diǎn)總數(shù)和自由度數(shù)的乘積。如平面問題，自由度為2，若有n個節(jié)點(diǎn)，則總剛度矩陣為:(2nx2n)階。若一個具有200個節(jié)點(diǎn)的小型平面問題，其總剛度矩陣的元素就有400x400=160000個，為一般中小型電子計(jì)算機(jī)的內(nèi)存容量所不允許的。因此，如何縮小總剛度矩陣所需的存貯單元，是有限元法程序編制中一個需要考慮的突出問題。我們可以根據(jù)總剛度矩陣的某些特性，尋求節(jié)省存貯量的途徑。1)總剛度矩陣是對稱陣。

利用對稱性，我們只需要存貯總剛度矩陣的上三角形部分（或下三角形部分）中的元素。

2)當(dāng)合理編排節(jié)點(diǎn)編碼時，總剛度矩陣可呈帶狀的稀疏陣。其中有不少元素為零，而非零元素對稱地分布于主對角線的兩旁，形成一帶狀陣。下面介紹兩種常見的壓縮存貯方法：(1)半寬帶存放：仍以圖2-16為例，分析與該圖相應(yīng)的總剛度矩陣，其元素的分布表示于圖2-18中

當(dāng)把對角線向右平移到第十個元素時，剩下的就全部是零元素了。

這是因?yàn)楣?jié)點(diǎn)編排時，就一個單元而言，節(jié)點(diǎn)編碼最大相差五個號碼，而每個節(jié)點(diǎn)又有兩個方向的位移所致。

我們把自對角線開始向右側(cè)平移至最后一根帶有非零元素的斜線為止，兩線之間一行內(nèi)包含的元素個數(shù)（包括兩線之間的零元素）稱為最大半帶寬。

而剛度矩陣中每一行的半帶寬均取決于一個單元中節(jié)點(diǎn)編號差和節(jié)點(diǎn)的自由度（即有幾個方向的位移）。

如以NBD表示最大半帶寬，以d來表示各單元中最大節(jié)點(diǎn)差值。（對于三角形單元，也就是相鄰節(jié)點(diǎn)編號的最大差值），則對具有兩個自由度的平面問題，最大半帶寬的計(jì)算公式為：NBD=（d+1）x2

對于以上所述的這樣一個對稱的，具有NBD半帶寬的總剛度矩陣，由于它在半帶寬以外有許多零元素，而對角線一側(cè)的元素又和另一側(cè)元素對稱相同，因此為了節(jié)省計(jì)算機(jī)內(nèi)存，可以采用長方形壓縮存儲法只將半帶寬內(nèi)的元素存放起來，如圖（2-18b）所示，這種存貯總剛度矩陣的方法也稱為“半帶寬存貯”或“等帶寬存貯”。圖2-18NEQNBDNBD

(a)

(b)

如節(jié)點(diǎn)總數(shù)為n，而每一節(jié)點(diǎn)有兩個方程，則方程總數(shù)為：NEQ=2n

我們可以用一個兩維數(shù)組SK[1：NEQ;1：NBD]即圖（2-18b）來存放總剛度矩陣中必須保存的元素。總剛度矩陣[K]和長方形壓縮存貯的數(shù)組[SK]之間有如下的對應(yīng)關(guān)系：矩陣[K]數(shù)組[SK]對角線第一列r行r行r行s列元素r行(s-r+1)列元素由[K]形成[SK]時，元素的行號相同。新的列號等于原先的列號減去行號加1，即:新列號=原列號-行號+1

為了節(jié)約內(nèi)存，我們力求減小半帶寬NBD，這就是要求在編排節(jié)點(diǎn)號碼時，應(yīng)使相鄰節(jié)點(diǎn)的節(jié)點(diǎn)差值d盡可能小。

[SK]右下角三角塊中的元素不和[K]中的任何元素對應(yīng)，那些元素應(yīng)是零。

考慮到總剛度矩陣[K]中各行的半帶寬并不相同，有時，由于結(jié)構(gòu)幾何形狀等原因，某些行的半帶寬特別大，而其它行又較小，這種情況下如以NBD為半帶寬（即等帶寬）存貯，就可能把許多零元素也存貯起來，這對節(jié)省計(jì)算機(jī)存儲量來說是不利的。(2)一維壓縮存貯法一維壓縮存貯法是將總剛度矩陣[K]的下三角形中每一行從第一個非零元素開始，逐行存放入一維數(shù)組K[1：S]中，S是元素的總個數(shù)）舉例如下，設(shè)有一對稱正定陣：在一維數(shù)組K[1：13]中依次存放為：

5，3，4，6，0，3，7，2，8，9，0，0，8

共13個數(shù)。

但是，僅使用這樣一個一維數(shù)組并不能將元素在[K]中的位置確定下來。為此，還須將主對角線上的元素在一維壓縮存貯中的序號用另一個一維數(shù)組N[1：2n0]存放起來（n0是節(jié)點(diǎn)總數(shù)）。

如對于上述矩陣數(shù)組N[1：6]中存放的是：1，3，6，8，9，13

即指出，在一維數(shù)組K[1：13]中，第1，第3，第6......個元素是對角元，顯然，這兩個數(shù)組完全確定了各元素在[K]中的位置。2.4方程求解1.引入位移約束條件

上一節(jié)，我們建立了有限元法的基本方程式：

有了這個方程還不能立即求解節(jié)點(diǎn)位移，因?yàn)榈浆F(xiàn)在為止，我們還沒有考慮到彈性體的幾何邊界條件，即邊界位移的約束條件.

很明顯，如彈性體的邊界沒有位移約束，則在外載荷的作用下，它將有產(chǎn)生剛體運(yùn)動的可能性，反映在基本方程上，其系數(shù)矩陣[K]將是一個奇異陣（對應(yīng)的行列式的值等于零）。逆矩陣不存在，方程將具有不定解。

在進(jìn)行應(yīng)力分析時，為了使解具有唯一性（即排除剛體運(yùn)動），必須根據(jù)彈性體具體的邊界位移約束條件，對基本方程加以處理，方能求解。1）引入約束條件的原因

（1）結(jié)構(gòu)實(shí)際上可能存在若干約束條件。它們應(yīng)該加以考慮，否則計(jì)算的結(jié)果將與實(shí)際不符。如一端固定，一端鉸支的靜不定梁，如圖2-19所示。圖2-19節(jié)點(diǎn)1，2，3及15既不允許有x方向的位移也不允許有y方向的位移。(2)有些約束是由于考慮到結(jié)構(gòu)和載荷的對稱性（或反對稱性）可取其中一部分作為計(jì)算對象而附加的。

如圖2-20所示一對角受壓的方形薄板。圖2-20

由于結(jié)構(gòu)和載荷的均對稱，可以只計(jì)算薄板的四分之一。

圖2-20b因?yàn)樽冃螌ΨQ于對角線，水平對角線nn‘上不可能有y方向的位移，所以附加的約束條件應(yīng)是：

（節(jié)點(diǎn)4，5，6為垂直的可動鉸支座）；

垂直對角線mm‘不可能有x方向的位移，所以附加的約束條件應(yīng)是:

（節(jié)點(diǎn)1，2，4為水平的可動鉸支座），薄板的中心點(diǎn)O，任何方向的位移均為零，故節(jié)點(diǎn)4處為固定鉸支座。2）引入約束條件的處理：

（1）零位移約束條件的處理

舉例說明，對于在剖分后有n個節(jié)點(diǎn)的彈性體，假定第n個節(jié)點(diǎn)處有約束，其位移為零，

即：un=0

vn=0

此外，還已知第（n-1）個節(jié)點(diǎn)沿y軸方向有約束，其y方向位移為零，

即：vn-1=0

（這對平面問題來說，是不產(chǎn)生剛體運(yùn)動的最低限度的邊界位移約束條件）。則應(yīng)對基本方程作這樣的處理：

把剛度矩陣[K]的最后三行和最后三列劃去，得到一個（2n-3）階的方陣，稱為總剛度矩陣的縮聚或降階。研究表明，降階以后的總剛度矩陣[K]是一個非奇異的對稱正定矩陣。

相應(yīng)地，分別將節(jié)點(diǎn)位移列陣{δ}和節(jié)點(diǎn)載荷列陣{R}的最后三行劃去，得到（2n-3）維列陣和

。

這樣，經(jīng)過零位移邊界條件處理后的基本方程就成為:

這樣得到的結(jié)果，不但滿足平衡條件和相容條件，而且滿足全部邊界條件，按上面的假定及處理方法，方程的具體形式是：式中Kij是第i行，j列的元素。

如果零位移的節(jié)點(diǎn)編號不在最后而在中間，也可以用同樣的方法劃去想相應(yīng)的行和列，而得到降階后的式（2-25）：

但是在計(jì)算機(jī)的計(jì)算程序中，我們采用的方法是使總剛度矩陣[K]中要劃去的行和列除對角線元素充成1以外，其余的元素均充成零。

如此得到的矩陣不降階，仍為2n階方陣，記為，并且也是一個非奇異的對稱正定矩陣。

與此同時，將節(jié)點(diǎn)載荷列陣中相應(yīng)的元素也充成零，如此得到的2n維列陣記為,因而基本方程就變?yōu)椋海?-26）

式中

為包含所有節(jié)點(diǎn)位移的列陣。顯然，在求解式（2-26）時，原來節(jié)點(diǎn)位移為零的仍保持為零。式（2-26）與式（2-25）實(shí)際上是等價的。

按前面同樣的假定，式（2-26）的具體形式應(yīng)該是：可以明顯看出，如果將上式乘開，則后三行的結(jié)果是：

也即：

vn-1=0,

un=0,

vn=0,

表示了邊界位移約束條件，由此可見，式（b）與式（a）等價。

有幾個約束條件，就應(yīng)重復(fù)幾次上面的修改過程。

處理邊界節(jié)點(diǎn)零位移條件，還可以采用其它方法，這里不一一列舉。（2）位移不為零的約束處理

如果某些節(jié)點(diǎn)的位移值是已知的，如：

我們可將基本方程作如下修改：

將第i行的對角元改為1

第i行的其它元素改為零

將右端項(xiàng)改為

式（2-27a）i行i列這樣修改后，第i個方程變?yōu)椋簽榱吮３挚倓偠染仃嚨膶ΨQ性，以便于以后求解，我們把[K]的第i列也改為零，這時，節(jié)點(diǎn)的載荷列陣要作如下修改：i行i列式（2-27b）

當(dāng)將（2-27a）式和（2-27b）式乘開時，兩組代數(shù)方程組完全相同。

下面舉一簡例說明:

已知δ3=u0，由（2-27a）式：（節(jié)點(diǎn)位移和節(jié)點(diǎn)載荷按列陣中的位置標(biāo)號）。由(2-27b)式:

方程（b）等同于方程(a)。

(2-27b)式保持系數(shù)矩陣對稱，目的是使求解方便。2.基本方程求解---高斯消去法

在把位移約束條件引入基本方程（2-23）式后，問題就歸結(jié)為線性代數(shù)方程的求解了。

一般來說，求解的方法可分成兩大類。

一類是直接法，本節(jié)要介紹的高斯消去法就是其中的一種；

另一類是迭代法，它是按照一定的迭代程序，由假設(shè)的初始值，通過多次迭代去逐步逼近方程的精確解。

直接法和迭代法相比，具有速度較快的優(yōu)點(diǎn)，但所需的存貯量比迭代法大，在計(jì)算機(jī)存貯量許可的條件下，通常都用直接法求解。

下面，我們通過一個簡單的算例，說明消去法的基本思想和步驟。求解線性代數(shù)方程組：

(a)

寫成矩陣形式為：

首先，可將（a）式第一式中X1的系數(shù)化為1，為此，用2去除第一式的兩邊,得到：

(b)

其次，我們把（b）式作為主導(dǎo)方程，使（a）的第二、三式X1的系數(shù)化為0

為此，只需用（+1）乘式（b），被第二式減，用（+3）乘式（b）被第三式減，如此，線性代數(shù)方程組可化為:寫成矩陣形式為：同理，可將（c）中第二式X2的系數(shù)化為1，得到:

(d)以（d）式為主導(dǎo)方程，將（c）式第三式X2的系數(shù)消去，得：最后得到下列形式的線性代數(shù)方程組：寫成矩陣形式：其中系數(shù)矩陣的特點(diǎn)是:

對角線系數(shù)為1，對角線下側(cè)的元素全部為零，這樣的方陣稱為

上三角陣。由（e）式第三式得：代入（e）第二式得：將x3及x2代入（e）第一式得：

至此，可把消去法求解線性代數(shù)方程組歸納成兩個主要步驟：

第一步是通過“化1消零”把線性代數(shù)方程組的系數(shù)矩陣化為上三角陣，這一過程稱為消元。

當(dāng)這一步驟完成時，方程組中最后一個未知數(shù)已求得；

第二步稱為反代，從最末第二個方程開始，將已求出的未知數(shù)代入，逐個求出上一個未知數(shù)。前面已指出：

總剛度矩陣[K]是一個對稱正定陣。其主對角線上所有元素均為正值。且有:

對于對稱正定陣，其逆矩陣一定存在。因此，以對稱正定陣[K]為系數(shù)的線性代數(shù)方程組：必有唯一解。（1）消元和反代過程：用矩陣表示：當(dāng)用高斯消元時，可將上式第一方程兩邊各除以K11，

使第一方程成為：（g）為主導(dǎo)方程。上式兩邊：

乘以K21后與（f）式第二方程相減，

乘以K31后與（f）式第三方程相減，

乘以K41后與（f）式第四方程相減，

便將方程化作如下形式：此時方程組（h）的系數(shù)矩陣的第一行第一列元素為1外，第一列的元素全部為零。用矩陣形式表示：

可以證明，第一循環(huán)后，自第二行和第二列開始的方程仍是一對稱正定陣。

因此，按同樣的方法，可使第二方程X2前的系數(shù)為1，而第三方程及其以下各方程在X2項(xiàng)的系數(shù)為零，依次類推，消元過程中的系數(shù)矩陣可表示如下：用矩陣形式表示消完元后的方程：1.，最后一個右端項(xiàng)的值就是最后一個未知量的解；2.從倒數(shù)第二個方程開始回代，每一個方程中都只有一個未知量，依次可求得全部節(jié)點(diǎn)位移值{}；（2）遞推公式：

上面的消元反代作法同樣適用于系數(shù)矩陣為n階的線性代數(shù)方程組，由上過程可知，消第一元時，系數(shù)矩陣從f1變到f2，也就是說，每消元一次，要作四件事（建立四個遞推公式）。（1）修改非主導(dǎo)方程的系數(shù)項(xiàng)

（2）修改主導(dǎo)方程的系數(shù)項(xiàng)

（3）修改非主導(dǎo)方程的右端項(xiàng)

（4）修改主導(dǎo)方程的右端項(xiàng)

其中,相應(yīng)于(2)(3)(4)項(xiàng)的遞推公式容易建立,唯有非主導(dǎo)方程系數(shù)項(xiàng)的處理涉及到有關(guān)元素如何從剛度矩陣[K]中過渡到[SK]中的問題。遞推公式：值得注意的是：1.每消元一次，需修改的方程只涉及到NBD個；

主導(dǎo)方程——1

非主導(dǎo)方程——2，3，4，...NBD2.隨著被修改方程的行數(shù)增加，方程中需修正的系數(shù)逐個減少（呈三角形，見下圖）。

3.在總剛[K]中的對角元，在[SK]中變?yōu)榈谝涣?。注?

（1）最后一個右端項(xiàng)就是最后一個未知量的解。

（2）反代從第NEQ-1個方程開始，倒推。反代結(jié)束。{R}陣是節(jié)點(diǎn)位移。

（3）每反代一個方程涉及該方程的NBD-1列元素和NBD個右端項(xiàng)。2.5應(yīng)力計(jì)算

全部位移求得后，各單元的節(jié)點(diǎn)位移確定后，單元內(nèi)的應(yīng)力可由公式（2-14）求得：其中應(yīng)力矩陣[S]=[D][B]

在計(jì)算程序中,計(jì)算某一單元的應(yīng)力時，首先要確定該單元的三個節(jié)點(diǎn)的六個位移分量是已求得的全部節(jié)點(diǎn)位移中的那幾個。相對應(yīng)的編號關(guān)系建立后，就可以把這些位移分量代入應(yīng)力分量關(guān)系式（2-14）求得σx，σy，τxy。

如仍以圖（2-16）為例，當(dāng)消元的反代過程結(jié)束后，節(jié)點(diǎn)位移分量全部求得：如欲確定單元（4）的應(yīng)力，按（2-14）式:需從全部位移分量中找出單元（4）相應(yīng)的位移分量。應(yīng)力分量求得后，主應(yīng)力和主方向按下式計(jì)算：2.6平面問題有限單元法小結(jié)

現(xiàn)在，我們來回顧一下用有限單元法按位移求解彈性平面問題的步驟和方法：（1）根據(jù)工程實(shí)際問題對計(jì)算精度的要求和計(jì)算機(jī)的容量，將平面彈性體剖分為一定數(shù)量的三角形單元，并對單元和節(jié)點(diǎn)進(jìn)行編號，選定右手坐標(biāo)系XOY（通常取水平X軸向右為正，垂直y軸向上為正），定出所有節(jié)點(diǎn)的坐標(biāo)值：

然后對每個單元按逆時針順序，確定i,j,m的具體節(jié)點(diǎn)編號，稱為單元節(jié)點(diǎn)信息。（2）根據(jù)節(jié)點(diǎn)坐標(biāo)值，用式（2-2），（2-3）標(biāo)出各個單元的幾何參數(shù)bi,ci,bj,cj,bm,cm和三角形面積Δ值。（3）再由這些數(shù)值，連同各單元的材料常數(shù)E和μ，按式（2-18b）或（2-20a）（平面應(yīng)力問題）或式（2-20b）（平面應(yīng)變問題）計(jì)算每個單元的剛度矩陣[K]e中的各元素。如為平面應(yīng)變問題，將式中的E換成

E/(1-μ2)和μ換成μ/(1-μ)。于是得到：

(4)按形成彈性體的總剛度矩陣[K]，即將每個單元剛度矩陣中的元素按節(jié)點(diǎn)編號的次序，填到總剛度矩陣相應(yīng)的行、列位置上去。如果在某個位置上填入若干個數(shù)，那就把它們迭加起來。（5）進(jìn)行載荷移置，算出彈性體總的節(jié)點(diǎn)載荷列陣{R}，從而得到以節(jié)點(diǎn)位移{δ}為未知數(shù)的基本方程式（2-23a）。（6）根據(jù)彈性體的幾何邊界條件，對上述基本方程進(jìn)行處理，得到式（2-25）（2-26）。（7）解上述線性代數(shù)方程組，求得全部節(jié)點(diǎn)位移。（8）用式（2-14）求出各單元的應(yīng)力分量σx，σy，τxy

，這些應(yīng)力在一個單元內(nèi)是常量，通常把它們當(dāng)作是作用在單元的形心上。另一種辦法是求出節(jié)點(diǎn)處的應(yīng)力；這時我們可認(rèn)為節(jié)點(diǎn)處的應(yīng)力等于環(huán)繞這個節(jié)點(diǎn)的若干個單元應(yīng)力的平均值。（9）根據(jù)單元形心或節(jié)點(diǎn)處的應(yīng)力分量，即可求得其主應(yīng)力及主方向。

在彈性體的不受外載荷作用的自由邊界上，一個主應(yīng)力的方向應(yīng)與邊界曲線相切，而另一主應(yīng)力的值應(yīng)為零。

當(dāng)彈性體在邊界上受到法向分布載荷作用時，一個主應(yīng)力方向同樣應(yīng)與邊界曲線相切，而沿法向的主應(yīng)力應(yīng)等于法向載荷的集度。

利用上述方法，也可以對用有限單元法求得的解答的精確度進(jìn)行估計(jì)。2.7矩形單元1.矩形單元的位移函數(shù)矩形單元有四個節(jié)點(diǎn)，按逆時針順序編號為1、2、3、4，設(shè)長為2a，寬為2b，厚度為t，為研究方便，使形心為坐標(biāo)原點(diǎn)，如圖所示。在平面問題中，每個節(jié)點(diǎn)有兩個位移分量，單元自由度為8，根據(jù)位移函數(shù)的完整性和收斂性要求，單元位移函數(shù)選為:在平行于x軸的直線上，位移分量是x的線性函數(shù)，在平行于y軸的直線上，位移分量是y的線性函數(shù)，故把上式稱為雙線性位移函數(shù)。在式(2-21)中代入節(jié)點(diǎn)位移和節(jié)點(diǎn)坐標(biāo)后，可解出待定系數(shù)。將這些數(shù)再代入式(2-21),可得:式中形函數(shù)為:則上圖的矩形單元變?yōu)樽鴺?biāo)系下邊長為2的正方形單元，節(jié)點(diǎn)1、2、3、4的坐標(biāo)依次為(-1，-1)、(1，-1)、(1，1)、(1，1)，如右圖所示。坐標(biāo)系稱為標(biāo)準(zhǔn)化坐標(biāo)系。式相應(yīng)變換為:令:于是位移函數(shù)可表示為:的變換可寫為:也可寫為:因?yàn)?通理可證:

坐標(biāo)變換式(2—27)與位移函數(shù)表達(dá)式(2—26)采用了相同的節(jié)點(diǎn)和相同的形函數(shù)，滿足這一特征的單元稱為等參元。由此可見，矩形單元是一種等參元。由平面問題幾何方程得:2．矩形單元的應(yīng)變矩陣、應(yīng)力矩陣與單元剛度矩陣其中:而:為3x8階矩陣,在xy

坐標(biāo)系中,表達(dá)式為:由于:即:于是,在坐標(biāo)系中:應(yīng)力矩陣也是3x8階矩陣。對于平面應(yīng)力問題，在坐標(biāo)系中，表達(dá)式為:單元剛度矩陣:

寫成分塊形式:而:所以:上式對應(yīng)的是平面應(yīng)力情況,對于平面應(yīng)變,請將式中的:E換成E/(1-μ2)μ換成μ/(1-μ)2.8ANSYS平面結(jié)構(gòu)計(jì)算示例

2.8.1問題描述如圖所示長方形板ABCD，板厚0.04m，孔半徑r=0.2m，E=210GPa，泊松比μ=0.3，約束條件：在長方形底邊AD約束全部自由度,BC邊施加垂直向下均布載荷g=