1、有限單元法及程序設(shè)計緒論1.力學(xué)分析方法:解析法,數(shù)值法有限元法一一實際結(jié)構(gòu)形狀和所受載荷比較復(fù)雜，大多用解析法很困難，因而數(shù)值法得到不斷發(fā)展，隨著電子計算機的進步,而發(fā)展起來的一種新興的數(shù)值分析方法.2基本步驟：（1）結(jié)構(gòu)離散化：將結(jié)構(gòu)從集合上用線或面劃分為有限個單元。（2）單元分析：導(dǎo)出單元的節(jié)點位移和結(jié)點力之間的關(guān)系（單元剛度矩陣）。（3）整體分析：將各單元組成的結(jié)構(gòu)整體進行分析，導(dǎo)出征個結(jié)構(gòu)點位移與結(jié)點力之間的關(guān)系。3程序設(shè)計的步驟：（1）提出問題，擬定解決方案（2）構(gòu)造數(shù)學(xué)模型（3）畫出程序流程圖（4）編寫程序（5）編譯調(diào)試程序（6）試算驗證程序4.根據(jù)國家標準（GBT526-89）

2、規(guī)定的程序流程圖標準化符號及規(guī)定：a）圖表示程序流程圖的起點和終點；端點框數(shù)據(jù)框處理框b）圖表示數(shù)據(jù)信息的輸入和輸出；C）圖表示數(shù)據(jù)進行系列運算之前要完成的數(shù)據(jù)預(yù)置:d）圖表示判斷條件：e）圖表示各種處理功能，如數(shù)學(xué)運算方式等；f）圖表示流程的路徑和指向。第一篇桿件結(jié)構(gòu)的有限單元法及程序設(shè)計第一章平面桿件單元的有限單元法第一節(jié)有限單元法的基本概念1 .基本思路：先分后合（先單元分析，再整體分析）2 .基本概念：整體號：節(jié)點端點號按自然數(shù)1,2,3,（在整體坐標系xOy下）局部號：每一個單元始末用i,j標記（在單元的局部坐標旗系下，方向與整體坐標系一致）。氐kjl3 .Fe=ke6*其中：ke=

3、單元剛度矩陣，各元素為剛度系數(shù)L品島6=單元桿段位移列陣F,=兇單元桿端力列陣KA=P(1-7)kuknki3K=k21k22k23整體剛度矩陣A=忸1仍仍T位移列k31k32k33陣P=MiM?M3節(jié)點載荷列陣3.有限元位移法分析連續(xù)梁需要考慮的問題（1）剛度集成法：將（1-3）K擴階,擴大的兀素為0,得到單兀貢獻矩陣單元：1kn023如kjj0單元：120kn如30kijkj將單元貢獻矩陣想登加,形成整體剛度矩陣K=K"+K叱kn品I0123kj10kJ+kiJkji3kJkJ（2）兩端支承條件的引入先不考慮約束條件,得到整體剛度矩陣后,將其主對角線元素kti改為1,第i行，第j

4、列其余元索改為0,對應(yīng)的轂荷元素也改為0.（3）非結(jié)點荷載的處理利用等效結(jié)點荷載進行分析：各結(jié)點（包括兩端結(jié)點）加約束，阻止結(jié)點轉(zhuǎn)動,其約束力矩分別為交于該結(jié)點的各相關(guān)單元的固端力矩之和，順時針為正.去掉附加約束（相當(dāng)在各結(jié)點施加外力荷載Pe，其大小與約束力矩相同，方向相反）將兩部分桿端方矩餐加起來.第二節(jié)局部坐標系中的單元剛度矩陣1.一般單元設(shè)單元的彈性模量、截面慣性矩、截面枳分別為A/、4桿長為1。單元的i、j端各有三個桿端力、X.和血（即軸力、剪力和彎矩）和與其相應(yīng)的三個桿端位移不、v,萬，如圖17所示。圖中天oy為單元局部坐標系，取i點位于坐標原點，京軸與桿軸重合，規(guī)定由i至燈為支軸的

5、正方向，由支軸順時針旋轉(zhuǎn)90為9軸正方向。力和位移的正方向如圖1-7所示。1在此單元中，單元桿端力列陣和桿端位移列陣分別為F=x!YiMiXj%mJS=UiMaUjVj京?單元桿端力列陣桿端位移列陣為了導(dǎo)出一般單元桿端力與桿端位移之間的關(guān)系，我們分別考慮以下兩種情況。首先分析兩個桿端軸力K、欠】與軸向位移«、q的關(guān)系。根據(jù)胡克定律，有-EA_、EA_EA_X=(u.-u.)=ux-u(a)、151jXj=-(Ux-u)=-ux+UjJIv1Jz11IJ其次考慮桿端方矩M、Kij與桿端剪力反、*與桿端轉(zhuǎn)角同、和橫向位移又、匕的關(guān)系。根據(jù)結(jié)構(gòu)力學(xué)位移法的轉(zhuǎn)角位移方程，并按照本節(jié)規(guī)定的符號

6、和正負號，可得M=黑+生空1看+里司1I2111I2J1J6EITY_"EIX-下-2EI-6EI_4EI-v+q-vi+e、111I2J1J6EI萬12EI_JEI5YF力F(b)12EI13_6EI_12EI_6EI-1j21/J2J將(a)、(b)兩式合在一起，并寫成矩陣形式如下上式可簡寫成其中單元剛度矩陣為ke=EAT0010012EI6EI112EI6EIi36EI-PI24EITI36EIT-2EIEAEAEATEAT12EIY6EIF6EI下2EIT12EI丁6EI-F(1-18)EA6EI4EI"T"12EII36EIT6EIT-4EIT"

7、;EAT12EI-p-6EI6EIF2EI12EIT-6EIT-12EIT-6EI6EIT2EIT"6EIT-4EIr(1-19)2 .單元剛度矩陣的性質(zhì)（r、s取1至6）的物理意義是第s個桿端位移分每個元素代表單位桿端位移引起的桿端力，任一元素匕量等于1時，所引起的第r各桿端力分量值.是對稱矩陣,其元素&二h（r工s）.是奇異矩陣，它的元素行列式等于零,即|k|=0.具有分快性質(zhì).3 .軸力單元：只考慮軸向桿端位移和桿端力的單元第三節(jié)單元剛度矩陣的坐標變換上述單元剛度方程和單元剛度矩陣實在局部坐標系天oy中建立起來的，對于一般桿件結(jié)構(gòu)，分析時所劃分的各單元的局部坐標系顯然不

8、同。因此在研究結(jié)構(gòu)平衡條件和變形連續(xù)條件時，必須選定一個統(tǒng)一的坐標系xOy,稱為整體坐標系。同時，還必須把在局部坐標系中建立的單元剛度矩陣轉(zhuǎn)換為整體坐標系下的單元剛度矩陣。圖1-8a）、圖l-8b）分別表示單元©在局部坐標系又oy和整體坐標系xOy種的桿端力分量。為了導(dǎo)出整體坐標系中桿端力兀、Y-和局部坐標系中尺、1、%之間的關(guān)系，將X、匕分別向x.y軸上投影，可得a)%=Xcos2+丫sina1X=-Xxsina+Xcosaj式中，a表示由x軸到京軸之間的夾角，以順時針為正,M=Mb)同理，對于單元上端的桿端力可得Xj=Xjcosa+Yjsina*=-Xjsina+Ycosa>

9、;Mx=MiJJ-c)將a）、b）、c）合起來，并用矩陣形式表示，可得cosasiiia0000X-sinacosa0000X001000M,000cosasilla0Xj000-sinacosa0¥000001Mj.此式即為兩種坐標系中單元桿端力的變換式，亦可簡叮為(1-24)F*=TeFe(1-25)式中：Fe=（XYMx5XjYjMjf局部坐標系中的單元桿端力列陣Fe=xYM：斗MmF整體坐標系中的單元桿端力列陣cosasina0000-smacosa0000001000000cosasilla0000-sinacosa0000001T為正交矩陣,其逆矩陣等于其轉(zhuǎn)置矩陣單元坐標

10、變換矩陣（L26）-1T即Te=Te第四節(jié)單元未知量編碼為了便于編程計算，需要按一定規(guī)律對結(jié)點的位移分量編號。結(jié)構(gòu)的節(jié)點位移有自由結(jié)點位移和支座結(jié)點位移（亦稱支座結(jié)點位移）之分。自由結(jié)點位移是未知量。建立結(jié)構(gòu)整體結(jié)構(gòu)方程求解未知節(jié)點位移的方式有兩種：“前處理法”和“后處理法工xb)圖1-9用后處理法分析改結(jié)構(gòu)時，設(shè)所有點位移都是未知量，則結(jié)點位移列陣為（參看圖1-9）=Aa3a3a41r=U1V6JU2VM；113V3%：114VPrPiy、p田分別代表作用在結(jié)點i（i=1，2,3,4）上的水平力、豎向力和力偶。規(guī)定，結(jié)點力Pnc、Pxy的正方向與整體坐標系X、y的正方向相同，Pw以順時針指向

11、為正；結(jié)點位移的正方向與結(jié)點力的正方向一致。在求出各單元剛度方程之后，根據(jù)結(jié)點平衡條件和位移連續(xù)條件，可建立整個結(jié)構(gòu)的位移法方程PiP2P3"P4k；kj+K：K；其中14(1-38)(1-37)%001(1-39)畤+匕耳00耳k-j+Ki氐00耳耳為結(jié)構(gòu)的整體剛度矩陣，或稱為結(jié)構(gòu)的原始剛度矩陣。在建立整體剛改方程式（1-38）時,假定所有點位移都是未知量，相當(dāng)整體結(jié)構(gòu)無知座因而在外力作用下，除了彈性變形外，還有可能發(fā)生剛體位移，此時，各結(jié)點位移不能唯一確定。這說明式(1-39)為奇異矩陣，不能求逆，故利用(1-38)不能求結(jié)點位移。實際在圖2-9a)所示的剛架中，結(jié)點1、4為固定

12、點，因此結(jié)點位移是己知的，支承條件全為零。將該支承條件引入到整體剛度方程，得耳p2P3P4(1-40)可以分為兩組方程，一組是它可以求結(jié)未知結(jié)點位移2、A>o令一組是稱為反力方程。利用式(1-41)求出結(jié)點位移2、Ag并代入上式后,便可計算未知的支座反力。對于i般桿件結(jié)構(gòu)，都可以按上述步驟進行分析。無論結(jié)構(gòu)有多少個結(jié)點位移分量，經(jīng)過調(diào)整其排列順序，總可以將它分為兩組：一組包括所有的未知結(jié)點的位移分量，以表示：另一組為支座結(jié)點位移分量，以工表示。相應(yīng)的，將全部結(jié)點分為兩組，與，相應(yīng)者為已知的結(jié)點力列陣，以PF表示；與K相應(yīng)者為支座結(jié)點力列陣，以仁表示。于是有fPf0=，此二r_Pr_與以上

13、分析方法相配合，將整體剛度矩陣K0中的各元素重新排列，則KOAgPO可寫成(1-43)展開上式得KfrAr=PfKppAp+KrfAf+KRRAR=PRKTrXVx.KR式(1-44)為“修正的整體剛度方程”，它與式(1-38)的區(qū)別在于引進了支承條件。后處理法：由單元剛度矩陣形成整體剛度矩陣,建立剛度方程后在引入支承條件,進而求解結(jié)點的位置位移的方法.前處理法：僅對未知的自由結(jié)點位移分量編號，得到的結(jié)點位移列陣中不包含已知的約束結(jié)點位移分量.圖1-10所示的具有組合結(jié)點的剛架劃分為三個單元，其編號為、，各桿之間的笳頭表示局部坐標系的正方向，剛架結(jié)構(gòu)編號為15。下面考慮各單元結(jié)點位移分量編號。

14、采用“先處理法”作如下規(guī)定：僅對獨立的位移分量按自然數(shù)編號，稱為位移號。若某些位移分量由于聯(lián)結(jié)條件或直接桿軸向剛性條件的限制彼此相等，則編號相同。在支座處，由于剛性約束而使某些位移分量為零時，此位移分量編號為零。因此圖因此編號如下桿數(shù)單元編號單元結(jié)點編號單元位移分量編號始端末端AB12000123BC23123456DC54000457在計算程序中，單元兩端結(jié)點號可采用二維數(shù)組JE(i,c)表示，稱為“單元兩端結(jié)點號數(shù)組”。JE(l,c)=單元始端的結(jié)點號7£(2戶)=單元末端的結(jié)點號在本例中JE(1,1)=1,JE(2,1)=2JE(1,2)=2,JE(2,2)=3JE(1,3)=

15、5,JE(2,3)=4任意結(jié)點位移分量的位移號可用二維數(shù)組JN(iJ)表示，稱為“結(jié)點位移號數(shù)組”JN(1,j)=結(jié)點j沿x方向的位移號JN(2,j)=結(jié)點j沿y方向的位移號JN(2,j)=結(jié)點j角位移的位移號在本例中，對第三結(jié)點而言，JN(1,3)=4JN(2,3)=5JN(3,3)=6將單元©始端及末端得位移號排成行(始端在前)，此數(shù)碼為“單元定位數(shù)組”，利用它可方便的形成桿端位移與相應(yīng)結(jié)點位移間的協(xié)調(diào)條件。它的展開式為mc=(叫,叱,皿)。式中，d個元素時叫分別是單元的兩端位移分量所對應(yīng)的位移號數(shù)值。在本例中,in1=(0nr=(1m3=(00043)6)7)第五節(jié)平面結(jié)構(gòu)的整

16、體剛度矩陣在進行了單元分析得出單元剛度矩陣之后，需要進行整體分析。以“先處理法”為例,將離散單元.重組合成原結(jié)構(gòu)，使其滿足結(jié)構(gòu)結(jié)點的位移連續(xù)條件和力的平衡條件，從而得到修正的結(jié)構(gòu)剛度方程，即前面給出的式（144）KffZf=Pf式中：K”稱為修正的結(jié)構(gòu)整體剛度矩陣：F、Pf分別為Hrti結(jié)點位移與自由結(jié)點荷載列陣。當(dāng)已知計算對象為自由結(jié)點位移分量而不至引起誤解時，式（144）也常稱為整體剛度方程，Kff簡稱為整體剛度矩陣,f、Pf分別簡稱為結(jié)點位移列陣與荷載列陣。為了書寫方便，下表常常略去。有單元剛度矩陣集成整體剛度矩陣，通常采用“直接剛度法”，把計算步驟分為兩步，首先求出各單元貢獻矩陣，然后

17、將它們疊加起來，得出整體剛度矩陣。然而在實際電算中，不便采用，原因是在計算中需要先將所有單元的貢獻矩陣K，都保存起來，K，的階數(shù)與整體剛度矩陣K相同，這就占用了大量存儲容量，因此實際運算中，采用“邊定位，便累加”的方法。其原理沒有變，而且結(jié)果相同。2(1二,3?-3(4空6).*飛©尸由7)八八1(0.0,0)A口5(0,0,0)圖1-10圖140中單元的單元矩陣為000123心fK12E5fK22心TK32145fK23式中單元剛度矩陣的上面和右側(cè)標記了單元結(jié)點位移分量編號.因為整體剛度矩陣個元素是按位移分量編號排列的，按先處理法，單元剛度矩陣中對應(yīng)于分量編號為零的元素不進入整體剛

18、度矩陣，非零編號指明了其余各元素在整體剛度矩正中的行、列號。所以A一4單元的剛度矩陣為123456各單元在K中的位置為單元的剛度知陣為K44曦一K75按以上定位方法，將三個單元的剛度矩陣有關(guān)元素一道整體剛度矩陣對應(yīng)位置，得到k%+k；心+弟扃上+匕心+舄心+匕K匕舄匕匕匕00在世機電算時,采用將K置零，這是K=07x7；將k中的相關(guān)元素，按照“對號入座”，累加到K；將k中的相關(guān)元素，繼續(xù)按照“對號入座”，累加到K；將k中的相關(guān)元素，繼續(xù)按照“對號入座”，累加到K,整體剛度矩陣最后完成。例13求圖1-12所示剛架的整體剛度矩陣Ko設(shè)各桿截面尺寸相同。A=0.5nr1=1/24m2£=3

19、X108解（1） 整理數(shù)據(jù)并進行編號。=300xl041=100x1041 16EI0八.八412EA.八4=30xl04=12xl04I2I3（2） 求局部坐標系中單元剛度矩陣1?。由于單元、的尺寸完全相同，故有!？=!?,可直接利用公式(1-19)求得3000001230_030100k1=k2=-300000-12-3003050-30000300000-12-30012-30030500-30100（3）求整體坐標系中單元剛度矩陣k。按式（1-32）和式（卜33）,求得各單元在整體坐標系中的單元剛度矩陣，并將單元結(jié)點位移分量編號標記于上面與左側(cè)。3k1=0(3)112-30-12-30

20、300-3001300-3003-3010030501230-12300-120-3030123030100-3050-3003000-3003005030100-12-3012-30xlO43050-30100IO4形成整體剛度矩陣Ko采用本節(jié)中介紹的方法建立結(jié)構(gòu)整體剛度矩陣如下3120-3003123030xlO4-30302005003050100圖1-12第六節(jié)非結(jié)點荷載處理為分析平面結(jié)構(gòu)而建立的整體剛度方程，反映了結(jié)構(gòu)的結(jié)點荷載與結(jié)點位移之間的關(guān)系。作用在結(jié)構(gòu)上的荷載除r直接作用在結(jié)點上的荷載匕之外，還有作用在桿件上的分布荷載等。這些非結(jié)點荷載應(yīng)轉(zhuǎn)換成等效結(jié)點荷載匕。將匕和、疊加，課的綜合結(jié)點荷載(總結(jié)點荷載)p“其下標c通常可略去不寫，即P=Pa+Pe(1-46)直接作用在結(jié)點上的荷載,可按其作用方位直接加入P之中.而等效結(jié)點荷載的計算步驟如下。第一步：在局部坐標系下，求單元的固端力F/。F/=|xflYfi取羽片MjT（1-47）FqLJ