1、第二章結(jié)構(gòu)矩陣分析由于有限元方法起源于力學(xué)中的結(jié)構(gòu)分析，本章的作用是通過三個典型問題說明有限元方法應(yīng)用于結(jié)構(gòu)分析時的一般步驟，并借此了解有限元方法的一些基本概念。2-1平面桁架（直接法，結(jié)構(gòu)矩陣分析中常用的力法，處理靜定問題，位移法，可處理靜定&靜不定）本節(jié)討論的對彖是圖2-1所示的平面桁架。組成桁架的各桿為等截面直桿，外載荷p直接作用于桿的較接點(diǎn)（結(jié)點(diǎn)）上。為簡單起見不妨設(shè)各桿的截面積均為A,材料的彈性模量均為E。我們可按下述步驟求得桁架的變形和內(nèi)力。22圖2-1圖22圖2-1圖2-21、結(jié)構(gòu)的離散化對結(jié)點(diǎn)及單元編號取組成桁架的每根桿為一個單元（該問題本身為一離散結(jié)構(gòu)的力學(xué)問題），以，加以編

2、號；取桿的餃接點(diǎn)為結(jié)點(diǎn)，以1、2、3加以編號（總體結(jié)點(diǎn)序號）。如圖2-2所示，即：我們所討論的桁架包扌舌三個單元、三個結(jié)點(diǎn)。各單元（桿）僅在結(jié)點(diǎn)處連接。2、建立總體坐標(biāo)系并確定結(jié)點(diǎn)坐標(biāo)和自由度為了描述結(jié)構(gòu)的平衡需要建立一個坐標(biāo)系，稱為總體坐標(biāo)系，以區(qū)別于以后出現(xiàn)的“局部坐標(biāo)系化總體坐標(biāo)系的選擇原則上不受限制，但希望使用方便。本節(jié)所選的總體坐標(biāo)系示于圖2-2,坐標(biāo)原點(diǎn)與結(jié)點(diǎn)1重合。以U,V分別表示沿X,y方向的位移分量，p,q分別表示力沿x,y軸的力分量（投影）。在總體坐標(biāo)系中各結(jié)點(diǎn)的坐標(biāo)為：y2）=（a,a）（x3yj它們將作為程序的輸入數(shù)據(jù)（幾何參數(shù)）。每個結(jié)點(diǎn)有兩個自由度，對結(jié)點(diǎn)1、2、3

3、分別為圖2-3ill.T、U2,V2T、I：V3T若暫時不考慮支承約束條件，整個結(jié)構(gòu)的結(jié)點(diǎn)自由度為illViU2v2U3V3T3、單元分析（建立結(jié)點(diǎn)力與結(jié)點(diǎn)位移之間的關(guān)系）取一個一般性的單元，設(shè)它的兩個結(jié)點(diǎn)在結(jié)構(gòu)中的編號為1,J（單元內(nèi)部的結(jié)點(diǎn)序號）。由材料力學(xué)可知，桿的軸向剛度為EA/L。其中L為桿的長度：L=J（x廠就+（為一）2（1）單元局部坐標(biāo)系現(xiàn)選取一典型單元對其進(jìn)行單元分析，對所分析的單元按如下方式建立一個坐標(biāo)系：原點(diǎn)：與結(jié)點(diǎn)1重合，x軸：沿i,j方向，y軸：與x軸垂直。如圖2-3所示。這個坐標(biāo)系只屬于一個單元，故稱為單元局部坐標(biāo)系，不同單元的單元局部坐標(biāo)系一般是不相同的。在單元局

4、部坐標(biāo)系中可以規(guī)定：結(jié)點(diǎn)自由度Ujv；T,UjVjT；單元結(jié)點(diǎn)自由度衛(wèi)=uJUjVjT。（2）局部坐標(biāo)系中的單元剛度矩陣在外載荷作用下，結(jié)構(gòu)發(fā)生變形，單元必受到來自結(jié)點(diǎn)的作用力。桁架中的桿只承受軸向力S,大小與桿的軸向伸長AL成正比EALEAL在局部坐標(biāo)系中這種特性可以得到清楚的表述（這一點(diǎn)也是引入局部坐標(biāo)系的理由之一）。若以pi,qi,pj,q3分別表示結(jié)點(diǎn)i,j作用于單元的力在x,y軸上的投影，由號單元的靜力平衡有（圖23）有，EA,PjUj(2亠(2亠1)EA,Pi=-s=-y-(i*_uj)q；=q；=o用矩陣的形式町以寫成若引入單元廣義力矢量：=p若引入單元廣義力矢量：=p：q；p；

5、q；T則上式可縮寫為0-10|q0000p：L一】010000Uxvj0(2-1-2)0000其中10-1000Ujcosasina其中10-1000Ujcosasina-sinacosacosasina-sinacosaIvj(2-1-3)稱為局部坐標(biāo)系中的單元剛度矩陣，它只與桿的幾個參數(shù)E、A.L有關(guān)，與桿的方位無關(guān)）（3）坐標(biāo)變換局部坐標(biāo)系中的單元剛度矩陣公式簡捷。但不同單元的局部坐標(biāo)系一般不同，為了研究結(jié)構(gòu)整體的平衡，必須將結(jié)點(diǎn)給單元的力以及相應(yīng)的單元剛度矩陣轉(zhuǎn)換到統(tǒng)一的坐標(biāo)系總體坐標(biāo)系。在總體坐標(biāo)系中單元結(jié)點(diǎn)自由度衛(wèi)=*ViUjVjT結(jié)點(diǎn)給單元的力r）=Pjq】pjQjT在圖2-3中

6、，x軸與x軸的夾角為a%杏y；-y：cosa=,suia=LL結(jié)點(diǎn)的位移分量的坐標(biāo)變換為00000000單元的位移分量的坐標(biāo)變換為Uxcosasina00U1UxV；-sinacosa00V,Vf001UjCOS6ZsinaUjUjJMJ00-siiiacosa%或縮寫為(2-1-4)(2-1-5)類似，/與r之間的轉(zhuǎn)換關(guān)系為(2-1-5)=t加由于t=cosa-siiicrsinacoscrt=cosa-siiicrsinacoscr(2-1-6)00000000是正交矩陣，因此(2-1-7)(2-1-7)也是正交矩陣。所以有T?=TF將(21二)、(2-1-5)代入(2-1-2)有Mr=k

7、lTfe從上式可得到fc=TTkrlTu=k(2-1-8)其中k=TTklT(2-1-9)稱為單元在總體坐標(biāo)系中的單元剛度矩陣。以后將會看到，(2-1-9)是一個具有普遍意義的公式。它表明，當(dāng)單元的自由度由一種形式換成另一種形式時，單元剛度矩陣只需進(jìn)行一次相似變換。對于平面桁架單元,將(2-1-3)、(2-1-6).(2-1-7)代入(2-1-9)可得到更便于應(yīng)用的單元剛度矩陣公式cos2ak卡cosasina-cos2acosasinasin2a-cos2a一cosasina一cosasina-sin2a一cosasina一cosasina-sin2acos2ak卡cosasina-cos2

8、acosasinasin2a-cos2a一cosasina一cosasina-sin2a一cosasina一cosasina-sin2acos*acosQsinacosasinasina(4)具體結(jié)果由(2-1-10)可求得各單元的剛度矩陣的具體形式如下：45單元剛度矩陣為illVlU2V21,a=4i血-血4444近-近-逅4444-近-VI忑血4444-近-VI忑近_4444叱uJT.a=_0000k2_EA010-1a0000單元：單元自由度單元：單元自由度0-101-90單元剛度矩陣為4=空a單元：單元自由度仙VujVjT(2-1-10)(2-1-11)(2-1-12)10-10k“E

9、A0a-1(2-1-13)k“EA0a-1(2-1-13)請注意，單元剛度矩陣與單元自由度中位移分量的排列次序有關(guān)。如果改動這種排列次序,例如對號單元，將單元自由度次序由U1VU2V2F改為u?V2UiVT,必然導(dǎo)致剛度矩陣（2-1-11）元素位置的變動。（5）單元剛度矩陣的物理意義和特點(diǎn)設(shè)平面桁架單元在總體坐標(biāo)系中剛度矩陣的一般形式為由（2-1-8）.當(dāng)單元結(jié)點(diǎn)位移為1000T時，在單元各結(jié)點(diǎn)上施加的力剛好為單元剛度矩陣中的第一列：k：Xk：1k3：k（1對壓的其他各列也可做出類似的解釋。即單元剛度矩陣的每一列相當(dāng)于一組特定位移下的結(jié)點(diǎn)力，如表2-1所示。由圖24可以獲得更為直觀的理解。表2

10、-1平面桁架單元剛度矩陣的物理意義單元結(jié)點(diǎn)位移作用于單元的結(jié)點(diǎn)力1000Thik2ik3ikiiT0100Tkuk22k32ki2T0010T心k?3k33也T0001T也k?4k34kuT對結(jié)點(diǎn)對結(jié)點(diǎn)3；kuku圖2-4對圖2-4中的各種情況，據(jù)平面力系的平衡條件應(yīng)有kis+k3s=k2s+k4s=（xJ-xiK-（yJ-y1K=0“（s=14）這三個關(guān)系說明，k的四個行向量中只有一個線性獨(dú)立（四個元素有三個約束方程）。從以上分析可以看出，一般的單元剛度矩陣均具備以卞兩個特征。（對平面桁架單元而言，從（2-1-10）也可以得出這些結(jié)論）（i）單元剛度矩陣是對稱矩陣，這是線性系統(tǒng)互易定理的具體

11、體現(xiàn)。由于對稱性，對行向量或列向量兩者之一得到的結(jié)論，對另一個也適用。（ii）單元剛度矩陣是奇異矩陣。它的行向量（或列向量）線性相關(guān)，具有零特征值，detk=0o對平面桁架的單元剛度矩陣而言，它的四個行向量（或列向量）中只有一個線性獨(dú)立，而k有三個零特征值。這三個零特征值對應(yīng)的特征向量相當(dāng)于三種獨(dú)立的剛體位移模式：兩個平移，一個旋轉(zhuǎn)。這是我們在單元分析中不考慮位移約束條件的自然結(jié)果。4、總體剛度矩陣的組裝總體平衡方程將圖2-1所示的桁架中的支承約束以約束反力代替，如圖2-5所示。下面來建立平衡問題的有限元方程。（1）結(jié)點(diǎn)平衡條件作用于圖2-5每個結(jié)點(diǎn)上的外載荷、支座反力以及來自單元的力應(yīng)處于平

12、衡。對結(jié)點(diǎn)1:對結(jié)點(diǎn)2:y軸上的投影，則單元m給結(jié)點(diǎn)12R1YR3Y以示結(jié)點(diǎn)1作用于單元m的力在的力在x,y軸上的投影應(yīng)為一pj對結(jié)點(diǎn)1:對結(jié)點(diǎn)2:y軸上的投影，則單元m給結(jié)點(diǎn)12R1YR3Y眼-Pi(1)-P】=RiY-q-Qi3=oP-P2-卩:=-Q2(1)-q2=0-P3-P3(3)=0Rsy-Qs_q3=0圖2-5(2-1-(2-1- )（2-M9）用來“書寫“組裝總剛度矩陣的過程是簡單而明了的。但事實(shí)上不能這樣做,將所有町生補(bǔ)充零元素擴(kuò)充為匡L(fēng)極人地浪費(fèi)了寶貴的存貯空間，這些零元素僅起到使kL的元素在總剛度矩陣中就位的作用。實(shí)際上采用的是另一種方法。如呆已選定各結(jié)點(diǎn)位移在結(jié)構(gòu)總體自

13、由度中的排列次序為U1ViU2V2U3V3對每個單元在形成單元剛度矩陣的同時還形成了個定位數(shù)組LM,它將指出單元自由度的各分量在總體自由度中的序號，如下表：單元號單元自由度LM(1)LM(2)LM(3)LM(4)U1V1U2V2T1234112V2113V3丁3456UiViu3v3t1256單元剛度矩陣中第s行第t列的元素kJ川到總剛度矩陣的第LM（s）行LM（t）列即可。這一組裝總體剛度矩陣的方法被形象的稱為“對號入座”。從（2-1-20）可以看出總剛陣匡是奇異陣，它的六個行向量（或列向量）中只有三個線性獨(dú)立。這是尚未考慮位移約束條件，結(jié)構(gòu)的剛體位移未受到限制的必然結(jié)果。（4）引入位移約束

14、條件由圖2J,平面桁架的位移約束條件為U1=VU1=V=V3=0代入方程（2亠18）,得到(2-1-21)EAF1EAF1+14V24-VI4_V24a44運(yùn)444-100_000-忑-104近0Rx4000&Y運(yùn)nnU2sP4VuV20近1+10-15040100Ay,-101_它顯然可以分成兩個方程組返4V240EA4它顯然可以分成兩個方程組返4V240EA4+140、U2、pV2=U3.0(2-1-22)4EAy/24a/2-1U2Xx=V00kj06.求單元內(nèi)力桁架單元的內(nèi)力只有一個軸力S。由(2-1-1)S的大小和正負(fù)與P；相同。由(2-1-2)和(2-1-4)t=k，M=klTfe

15、取出它的第三行即得到EAr一cosaLL-siiiacosasina(2-1-24)EAr一cosaLL-siiiacosasina(2-1-24)對單元、求得內(nèi)力s于下表單元號單元結(jié)點(diǎn)位移內(nèi)力（以拉為正）小匚八PaPa丁0o（2J2+1）-TEAEA（2V2+1）Pa-Pa00TEAEA-p0000T0由于本節(jié)所討論的桁架是一個十分簡單的靜定桁架，用理論力學(xué)的知識即可得到各桿的內(nèi)力，結(jié)果相同。但是，若桁架為靜不定桁架且桿的數(shù)目有上千個，那么本節(jié)所討論的方法原則上不會遇到任何困難，我們要解決的課題將轉(zhuǎn)為如何管理有關(guān)的人量數(shù)據(jù)和如何解一個數(shù)下階的代數(shù)方程組這樣一些技術(shù)問題。2-3平面應(yīng)力問題常應(yīng)

16、變?nèi)切螆D2-8為一邊長為a、厚度為t的正方形薄板。其中AB邊固定，BC、CD邊自由，AD邊作用均布壓力q。對這一問題，有限元分析的步驟是：1、將ABCD劃分（離散）為8個三角形（單元），編號一。各單元僅在頂點(diǎn)（結(jié)點(diǎn)）較接，結(jié)點(diǎn)編號19。建立坐標(biāo)系后，不難定出各結(jié)點(diǎn)的坐標(biāo)（冷）。2、單元分析任取一個一般性的單元，如圖29所示。三個結(jié)點(diǎn)的編號為i,j,ko結(jié)點(diǎn)位移為（Ui,E）、（號號）、（uk,vfc）單元結(jié)點(diǎn)位移為也=hUj八vj圖u=+v=z4+cr5x+6Z6y（2-3-1）yyqj八vj圖u=+v=1%y廠=da2_1卷Yk.J可解出:&6為待定常數(shù)，在結(jié)點(diǎn)處應(yīng)有6Z-=da1-1其中

17、ax=ykl1Ay/、Uj=1%y廠=da2_1卷Yk.J可解出:&6為待定常數(shù)，在結(jié)點(diǎn)處應(yīng)有6Z-=da1-1其中ax=yklVkJ12Ay：yk1AYj=y3-ykiYk111y】旳yk=xky1-x1ykak=1Yiiy廣112A=detDA=ax+a;+akI】=右（6+g+cyh+2（aj+bjX+Cjyhj+右仇+btx+5*=N（x,y）i*+叫（兀y）uj+Nk（x,y）uk當(dāng)i?j,k的位置為逆時針排列時，2恒正，且等于三角形單元面積的兩倍。將這些結(jié)果代入(2-3-1)有類似可得到v=Nx(x,y)%+Nj(x,y)vj+Nk(x,y)vk可以合并成V1KNx00Nk000N

18、k0W0N0Nj0NkJVj0N0Nj0Nk_（2）單元的應(yīng)變、應(yīng)力利用（231）不難求得3LO慝NkOoNJJNooN*1No3LO慝NkOoNJJNooN*1No一一o0一內(nèi)舐Aax0內(nèi)-llrj一o0一內(nèi)0一舐Aax0內(nèi)2Ac】0C!qq0Cjobk0Cj0Cku=bu5Ckbk其中B(2-3-3)(24)在假定單元內(nèi)位移場u、v是x,y的一次函數(shù)的前提下，單元內(nèi)的應(yīng)變和應(yīng)力將是常數(shù),故這種單元又稱為常應(yīng)變?nèi)窃?。圖(b)(2-10)圖(b)(2-10)（3）為了在單元內(nèi)構(gòu)成均勻應(yīng)力場，必須在單元的各邊施加均布載荷，它們的合力一定作用在各邊的中點(diǎn)，如圖2-10(3)所示。再將各邊上的合力

19、平分到這邊的兩點(diǎn)結(jié)點(diǎn)，由圖240(b)不難得出Pi=片(Fix+f3x)=|(yj-y：R+(-%)rxy-(yk一兀-=$(yj-ykbx+(耳一Xjkxy類似可求得q】、Pj、qjpkqk并可合并寫成Pi10cqPi10cq：Cxb,*t=50CJq20CjbJpk%0CkQt.0ck=tABT=tABTEjBu根據(jù)單元剛度矩陣k的直觀意義，常應(yīng)變?nèi)窃膯卧獎偠染仃嚰礊?2-3-5)(4)單元和的邊上作用著均布的外載荷，可以把它們的合力平分到兩結(jié)點(diǎn)，如圖2-11所示。f2=：q&-X4)辦/5=fl丘5辦/5=fl丘5（5）為了組裝總體剛度矩陣，每個單元還應(yīng)形成一個數(shù)組LMo元素LM（1

20、）LM（6）分別為2、巧、巧、業(yè)、在總體自由度中的序號。由于單元的結(jié)點(diǎn)位移和結(jié)點(diǎn)力都是用總體坐標(biāo)描述的，單元剛度矩陣不必再進(jìn)行坐標(biāo)變換。3、組裝總體剛度矩陣和載荷向量。為了實(shí)現(xiàn)位移約束條件U?=v7=us=vs=0=v9=o這些自由度對應(yīng)的行和列可以不必組裝。總體平衡方程為kQ=f（2-3-6）其中,非約束自由度位移為Vu2v2u3U4v4u5v5u6V6T相應(yīng)自由度的載荷向量為f=o-f200000-（f+fj00ooT解方程（2-3-6）可得到各非約束自由度的位移，再由（2-3-3）可求得各單元的應(yīng)力。2-4結(jié)束語1、本章所討論的三具體問題盡管力學(xué)背景不同，但分析步驟人體相同。區(qū)別僅在于結(jié)

21、點(diǎn)參數(shù)不同，單元分析公式不同。而組裝總體矩陣的方法、處理約束條件的方法以及代數(shù)方程組的求解方法完全通用。2、從建立有限元方程的方法看，本章屬于直接方法一類。這種方法只在簡單情況下有效，但是物理意義清晰。盡管今天有限元方法已發(fā)展成為一種抽彖的數(shù)學(xué)分析方法，本章提供的直觀解釋仍然常常被用來說明一些問題。3、對于第一、二兩個問題（平面桁架和平面框架）本章求的將是精確解（真實(shí)解）。而對第三個問題（平面應(yīng)力問題），本章求得的只是近似解（位移場是假設(shè)的），這樣就自然引出以下兩個問題：（1）所求的近似解與真實(shí)解的接近程度如何？或者說，怎樣才能接近得更好一些？（2）在進(jìn)行單元分析時，曾作了兩條關(guān)鍵性的假設(shè)：（

22、i）單元內(nèi)位移u、v是坐標(biāo)x,y的一次函數(shù)：（11）各單元僅在結(jié)點(diǎn)處較接。這樣的假定有沒有依據(jù)？是否必需？這些問題將在后續(xù)的章節(jié)中討論。變分原理簡介1、引言變分涉及求以函數(shù)為變量的函數(shù)一一泛函的極值的問題。我們首先以三個非常簡單的極值問題為例，引出泛函和變分的概念。這三個問題是：最短路徑問題；最速下落問題：最小旋轉(zhuǎn)曲面問題。（a）最短路徑問題求連接給定兩點(diǎn)的最短平面曲線。坐標(biāo)系取為直角坐標(biāo)系，其中的一個點(diǎn)位于原點(diǎn)，另一個點(diǎn)的坐標(biāo)為（旳，力）。如果V-是一條通過（0,0）和（xi,刃）的曲線，I是兩點(diǎn)間曲線的長度，顯然XiI=J77?分耳2護(hù)(1)jV1+dx,(1)所求問題為，求函數(shù)y,使上述

23、積分值為最小。（b）最速下落問題求鉛直面內(nèi)一曲線，當(dāng)粒子在重力作用下從曲線的一點(diǎn)由靜止?fàn)顟B(tài)滑落到曲線上另一點(diǎn)時所用時間最少。建立如下圖所示的坐標(biāo)，起點(diǎn)位于原點(diǎn)，終點(diǎn)坐標(biāo)為yyds（xbyi）o任一時刻粒子運(yùn)動的速度航為勺莎。粒子從原點(diǎn)卞落到（xi,y】）所需時間為所求問題為，求函數(shù)y,使上述積分值為最小。（O最小旋轉(zhuǎn)曲面問題給定兩點(diǎn)和一條直線，點(diǎn)和直線在一個平面內(nèi)，求連接兩點(diǎn)的一條曲線，使曲線繞直線旋轉(zhuǎn)生成的曲面面積為最小。取直線為x軸，其中的一個點(diǎn)在y軸上，曲線端點(diǎn)的坐標(biāo)為（0,y0）和（xi,y】）。用S表示旋轉(zhuǎn)曲面的面積，則S=yds=o記所求問題為，求函數(shù)y,使上述積分值為最小。這三個

24、問題的共同特點(diǎn)是，在給定條件下，求函數(shù)尸f(x)使積分J/(禺y,/)加為最小。我們先來求解第一個問題。最短路徑對應(yīng)的積分式可表示為I=J11+子2dx0設(shè)所求曲線方程為y=f(x),而尸F(xiàn)(x)是不同于f(x)的任意一條曲線的方程。兩者之差為叩)=F(x)一幾嘰位于F(x)和f(x)之間的曲線可表示為2/=+切)，其中a是不依賴于x的參數(shù)。如果a足夠小，則W=陶+切)位于尸f(x)的鄰域內(nèi)。如果曲線fVTWTa)2尸f(x)和尸F(xiàn)(x)是給定的，則積分僅是a的函數(shù)。a=0處枳分取最小值的條件為aa)=0,即廠(0)=0。F(a)=j71atf(妨陸0HH1+y+a?(x)P在a=o處r(r(

25、o)=Ylf(x)dx對上式進(jìn)行分部積分在端點(diǎn)處差值4(龍)=F)一3)=0。由于譏龍)是任意的，故上式積分為零的的條件求解該微分方程，得求解該微分方程，得/y,=C,71+滬y=K,y=Kx+L代入邊界條件（Q）anda】，小）,得y=x,即所求曲線為過給定點(diǎn)的直線。下面進(jìn)一步求解第二個問題，導(dǎo)出求解泛函極值更一般的表達(dá)式。JG（場y,以）dx以上三個問題的一般表達(dá)式為兀,參照問題一的求解過程，對給定的曲線J0（禺的獷）dxF（x）和f（x）,位于兩者之間的其它曲線方程為2/=+切）,積分.僅是a7（a）=Jy+叩（龍），yf+arjx）d.(a)=的函數(shù)。枳分呎對求a導(dǎo)數(shù)，得(a)=（禺V

26、+切（妨，W-ra（R）dR作)u=y+a?(x),砂=/+a?/(R)du(、dvda=作)才】a）=/U,砂）必=J瓠(”)+瓠)必(0)=對上式進(jìn)行分部積分，得廠（0廠（0）=J駕祕如+r（o）J陽-鶴如since乃(龍)vanisheswhenx=Xqandwhenx=Xi.Zz(0)mustbezeroindependentlyoftheformofrj(x)therefore該偏微分方程稱作拉格朗口方程。Forthecurveofquickestdescent0=亠71+/v.Art.1(b)(2).yxd(d(t=0,LagrangesEquationbecomesyFortheminimuynsurfaceLagrangesEquationbecomesyFortheminimuynsurface0=十2=VT+V5!v.Art.1(c)(3).LagrangesEquationbecomesVI+yVI+y/2-二、變分符號