平面桿件結(jié)構(gòu)的有限單元法-濟(jì)南大學(xué)土木建筑學(xué)院-土木建筑…

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

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

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

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

5、正方向，由支軸順時(shí)針旋轉(zhuǎn)90為9軸正方向。力和位移的正方向如圖1-7所示。1在此單元中，單元桿端力列陣和桿端位移列陣分別為F=x!YiMiXj%mJS=UiMaUjVj京?單元桿端力列陣桿端位移列陣為了導(dǎo)出一般單元桿端力與桿端位移之間的關(guān)系，我們分別考慮以下兩種情況。首先分析兩個(gè)桿端軸力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ī)定的符號(hào)

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

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

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

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

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

11、為正；結(jié)點(diǎn)位移的正方向與結(jié)點(diǎn)力的正方向一致。在求出各單元?jiǎng)偠确匠讨螅鶕?jù)結(jié)點(diǎn)平衡條件和位移連續(xù)條件，可建立整個(gè)結(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）時(shí),假定所有點(diǎn)位移都是未知量，相當(dāng)整體結(jié)構(gòu)無知座因而在外力作用下，除了彈性變形外，還有可能發(fā)生剛體位移，此時(shí)，各結(jié)點(diǎn)位移不能唯一確定。這說明式(1-39)為奇異矩陣，不能求逆，故利用(1-38)不能求結(jié)點(diǎn)位移。實(shí)際在圖2-9a)所示的剛架中，結(jié)點(diǎn)1、4為固定

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

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

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

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

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

17、將它們疊加起來，得出整體剛度矩陣。然而在實(shí)際電算中，不便采用，原因是在計(jì)算中需要先將所有單元的貢獻(xiàn)矩陣K，都保存起來，K，的階數(shù)與整體剛度矩陣K相同，這就占用了大量存儲(chǔ)容量，因此實(shí)際運(yùn)算中，采用“邊定位，便累加”的方法。其原理沒有變，而且結(jié)果相同。2(1二,3?-3(4空6).*飛©尸由7)八八1(0.0,0)A口5(0,0,0)圖1-10圖140中單元的單元矩陣為000123心fK12E5fK22心TK32145fK23式中單元?jiǎng)偠染仃嚨纳厦婧陀覀?cè)標(biāo)記了單元結(jié)點(diǎn)位移分量編號(hào).因?yàn)檎w剛度矩陣個(gè)元素是按位移分量編號(hào)排列的，按先處理法，單元?jiǎng)偠染仃囍袑?duì)應(yīng)于分量編號(hào)為零的元素不進(jìn)入整體剛

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

19、X108解（1） 整理數(shù)據(jù)并進(jìn)行編號(hào)。=300xl041=100x1041 16EI0八.八412EA.八4=30xl04=12xl04I2I3（2） 求局部坐標(biāo)系中單元?jiǎng)偠染仃??。由于單元、的尺寸完全相同，故有!？=!?,可直接利用公式(1-19)求得3000001230_030100k1=k2=-300000-12-3003050-30000300000-12-30012-30030500-30100（3）求整體坐標(biāo)系中單元?jiǎng)偠染仃噆。按式（1-32）和式（卜33）,求得各單元在整體坐標(biāo)系中的單元?jiǎng)偠染仃嚕卧Y(jié)點(diǎn)位移分量編號(hào)標(biāo)記于上面與左側(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é)點(diǎn)荷載處理為分析平面結(jié)構(gòu)而建立的整體剛度方程，反映了結(jié)構(gòu)的結(jié)點(diǎn)荷載與結(jié)點(diǎn)位移之間的關(guān)系。作用在結(jié)構(gòu)上的荷載除r直接作用在結(jié)點(diǎn)上的荷載匕之外，還有作用在桿件上的分布荷載等。這些非結(jié)點(diǎn)荷載應(yīng)轉(zhuǎn)換成等效結(jié)點(diǎn)荷載匕。將匕和、疊加，課的綜合結(jié)點(diǎn)荷載(總結(jié)點(diǎn)荷載)p“其下標(biāo)c通常可略去不寫，即P=Pa+Pe(1-46)直接作用在結(jié)點(diǎn)上的荷載,可按其作用方位直接加入P之中.而等效結(jié)點(diǎn)荷載的計(jì)算步驟如下。第一步：在局部坐標(biāo)系下，求單元的固端力F/。F/=|xflYfi取羽片MjT（1-47）FqLJ