平面桁架程序的計算原理及程序編制

結構有限元程序設計基本原理——平面桁架程序的計算原理及程序編制5-01矩陣位移法5-02矩陣位移法算例5-03最小總勢能原理的應用5-04矩陣位移法的求解步驟5-05結構計算簡圖的數(shù)據(jù)結構5-06位移未知數(shù)的確定5-07單根桿件的分析5-08結構總剛度矩陣的形成5-09桿件內力的計算5-10能量原理和矩陣位移法平面桁架程序的計算原理5-1矩陣位移法桁架是由離散桿件組成的構架結構，桿件的端點借助于無摩擦的鉸連接起來。桁架主要靠各桿中的軸向拉力或壓力來傳遞作用在桁架節(jié)點上的荷載，桿件的任何彎曲均忽略不計。用有限元分析桁架時，桁架中的每根桿件都是一個單元，稱為桿單元。它是一維單元。不同分類：（1）平面桁架/空間桁架（2）靜定桁架/超靜定桁架求解方法力法：以力未知數(shù)，必須預先滿足平衡條件，然后通過連續(xù)條件求解未知力；超靜定基的選取。位移法：以位移為未知數(shù)，各桿件的變形由相連接的節(jié)點位移確定（變形協(xié)調條件），通過各個節(jié)點的平衡方程求出未知位移，再由位移計算出各桿件的內力；各節(jié)點的平衡方程也可由最小總勢能原理推導。以平面桁架結構分析的程序設計為例，介紹結構分析和程序設計的方法。5-2矩陣位移法算例如圖所示平面桁架，各根桿的截面積F相等，材料的彈性模量E相同，在兩個單位力P=1的外荷載作用下，用位移法計算節(jié)點位移和各桿內力。平面桁架結構1234PP節(jié)點總數(shù)：NW=4可動節(jié)點數(shù)：NU=2位移未知數(shù)總數(shù)：NDISP=2*NU=4注意：節(jié)點自由度排序和節(jié)點平衡方程相對應。根據(jù)虎克定律，對于任意一根桿，有：考慮單根桿件在位移下產生的內力。式中為桿的伸長，為桿的長度，為桿的內力，稱為單根桿的剛度(單元剛度陣)。單元局部剛度陣用RD來表示單元剛度陣，于是胡克定律可表達如下由單根桿件的變形幾何關系可得

角是桿件軸線與方向的夾角，由正向逆時針向轉至桿軸的角度為正。進一步有內力和節(jié)點位移之間的關系對于每根桿件（以兩端節(jié)點編號A和B定出角）應用上述公式，有應用上述各桿內力和位移關系后，便可建立以位移為未知數(shù)的節(jié)點平衡方程。平面桁架結構節(jié)點平衡方程為

注意：節(jié)點自由度排序和節(jié)點平衡方程相對應。將各桿的內力用位移表示的方程代入上式，有未知數(shù)為節(jié)點位移寫成矩陣形式，即正則方程。其中結構剛度矩陣而外力向量為：采用高斯消元法求得位移為利用每根桿的內力-位移關系計算桿內力5-3最小總勢能原理的應用總勢能由兩部分組成結構的彈性應變能外力由于結構變位所產生的勢能最小總勢能原理與位移法都是以位移為未知數(shù)使變形狀態(tài)預先滿足連續(xù)條件。現(xiàn)對上述例題采用最小總勢能原理進行求解。對于整個桁架應變能是所有桿件應變能的疊加，即由于桁架結構的桿件只能承受拉壓力，所以單根桿件的應變能為把上述值代入應變能表達式，得到由上式可見，公式中只有位移的二次項，也就是說是位移的一個二次齊次函數(shù)，或者說是一個位移的二次型。（位移的正定二次型，應變能總是正的）對于該平面桁架，有桿號12345Laaasqrt(2)asqrt(2)a△Lv1u1-u2v2(u1+v1)/sqrt(2)(v2-u2)/sqrt(2)現(xiàn)在計算外力勢能。外力產生勢能的原因是當節(jié)點發(fā)生位移時，外力要作功。所作功的負值便是它們具有勢能的改變量，如果取未變形位置外力的勢能為零，有平面桁架結構1234PP將和相加，得到總勢能。由于是位移的一次函數(shù)，總勢能就成為位移的二次非齊次函數(shù)。根據(jù)最小總勢能原理，在所有可能的位移狀態(tài)中，真正發(fā)生的位移狀態(tài)使總勢能為最小。即函數(shù)對自變量的偏微商為零，即式中現(xiàn)對各位移變量分別取偏微商后，得注：值得指出的是剛度矩陣中的系數(shù)只與結構本身的幾何形態(tài)和約束條件有關，而與外載無關。5-4矩陣位移法的求解步驟結構計算簡圖->（節(jié)點、單元編號，建立一個統(tǒng)一的坐標系等）分析節(jié)點位移的力學特性->確定位移未知數(shù)。建立每根桿件兩端位移和內力的關系->（單根桿件的剛度矩陣）根據(jù)每根桿件上的上述關系建立結構可動節(jié)點的平衡方程->（結構總剛度矩陣）求解平衡方程，得到節(jié)點位移。根據(jù)求得的位移，利用每根桿件位移與內力關系計算各桿的內力。5-5結構計算簡圖的數(shù)據(jù)結構完整而確切描述一個平面桁架結構的數(shù)據(jù)有三個方面：（1）結構本體描述數(shù)據(jù)（NW,IESG,NU,X,Y,HL,HR）（2）性質數(shù)據(jù)（EF）（3）荷載數(shù)據(jù)（PX,PY）NW為節(jié)點總數(shù)IESG桿件總數(shù)NU可動節(jié)點總數(shù)X,Y節(jié)點坐標HL,HR每根桿件兩端節(jié)點編號EF性質數(shù)據(jù)PX,PY外載荷數(shù)據(jù)5-6位移未知數(shù)的確定對于平面桁架，每個節(jié)點有兩個自由度，把第個節(jié)點的水平位移、垂直位移分別記為，，這樣結構共有2·NU個位移。位移的方向與坐標軸相同為正，以這些位移作為未知數(shù)，并排列成一個列向量，稱為結構總位移向量。第i個節(jié)點的位移在總位移向量中占i0+1和i0+2位置，而i0=2*i-2。5-7單根桿件的分析根據(jù)虎克定律，對于任意一根桿，有：式中為桿的伸長，為桿的長度，為桿的內力，稱為單根桿的剛度（單元剛度陣）。一、桿件在局部坐標系中的剛度矩陣

注意⊿L=uB-uA，并且由平衡關系得到可把FA、FB排成列向量{N}，uA，uB，排成列向量{u}，系數(shù)陣排成矩陣RD，即設桿端分別得到平行于xoy坐標軸的位移uA，vA，uB，vB；則桿件的伸長量為二、位移的坐標轉換把uA，vA，uB，vB

寫成列向量{V}，系數(shù)排列成行向量{B}T，上式可以寫成如下形式其中把⊿L代入上面的公式，得到此時，已用全局坐標系中的位移表達出桿的內力N。由于最終得到的平衡方程都是相對于全局坐標系建立的，上面所計算出的內力與桿軸方向一致。假設在桿端作用平行于桿軸的力FA、FB，則由平衡方程FA和FB在x軸、y軸方向的分量分別為：把、、、排列成列向量，則有由此令得到全局坐標系下的單元剛度矩陣5-8結構總剛度矩陣的形成假設平衡方程形式為展開上式可以寫成式中nn=2*NU；P1，P2，…，Pnn為加在節(jié)點上的外力。5-9桿件內力的計算對任意一根桿件，如果桿的兩端位移已知。這些位移在總位移向量中i0+1，i0+2，j0+1，j0+2的位置上，表示如下進而用下面公式計算結構內力平面桁架計算程序的編制矩陣位移法求解一般步驟：（1）結構計算簡圖。（節(jié)點、單元編號，建立一個統(tǒng)一的坐標系等）（2）分析節(jié)點位移的力學特性，確定位移未知數(shù)。（3）建立每根桿件兩端位移和內力的關系。（單根桿件的剛度矩陣）（4）根據(jù)每根桿件上的上述關系建立結構可動節(jié)點的平衡方程。（結構總剛度矩陣）（5）求解平衡方程，得到節(jié)點位移。（6）根據(jù)求得的位移，利用每根桿件位移與內力關系計算各桿的內力。程序流程平面桁架計算程序的總功能表示{依據(jù)平面桁架及計算模型的有關數(shù)據(jù)，進行結構靜力計算，并輸出計算結果}根據(jù)平面桁架的特性進一步展開為：{總體數(shù)據(jù)結構的設計}{輸入描述平面桁架本體結構和性質的數(shù)據(jù)}{輸出結構圖形}{計算各桿件的單元剛度矩陣}{形成結構總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項總外力向量}{回代求解，求出總位移向量}{打印位移}{計算各桿內力}{輸出各桿內力}根據(jù)平面桁架的特性進一步展開為：{總體數(shù)據(jù)結構的設計}{輸入描述平面桁架本體結構和性質的數(shù)據(jù)}{輸出結構圖形}{計算各桿件的單元剛度矩陣}{形成結構總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項總外力向量}{回代求解，求出總位移向量}{打印位移}{計算各桿內力}{輸出各桿內力}{總體數(shù)據(jù)結構的設計}平面桁架程序的總體數(shù)據(jù)結構設計，包括全局量（標識符）的存儲安排與說明。一個平面桁架結構的數(shù)據(jù)有三個方面：結構本體描述數(shù)據(jù)（NW,IESG,NU,X,Y,HL,HR）性質數(shù)據(jù)（EF）荷載數(shù)據(jù)（PX,PY）NW為節(jié)點總數(shù)IESG桿件總數(shù)NU可動節(jié)點總數(shù)X,Y節(jié)點坐標HL,HR每根桿件兩端節(jié)點編號EF性質數(shù)據(jù)PX,PY外載荷數(shù)據(jù)MODULEPDATAREAL,PUBLIC::BAR_HL(100),BAR_HR(100),&BAR_TENSION(100),BAR_EF(100)REAL,PUBLIC::NODP_X(80),NODP_Y(80)ENDMODULE

REALB(4),DP(120),R(1200),C(4,4)INTEGERII(4),V(120),NW,IESG,NU{總體數(shù)據(jù)結構的設計}4‘NW4‘IESG2‘NU01‘X;Y11100014100‘HI;HL;EF12100231001310010‘PX;PY10節(jié)點總數(shù)：NW=4可動節(jié)點數(shù)：NU=2位移未知數(shù)總數(shù)：NDISP=2*NU=4平面桁架結構1342{輸入描述平面桁架本體結構和性質的數(shù)據(jù)}結構本體描述數(shù)據(jù)（NW,IESG,NU,X,Y,HL,HR）性質數(shù)據(jù)（EF）READ(1,*)NWREAD(1,*)IESGREAD(1,*)NUNDISP=NU+NUWRITE(2,*)'NW=',NW,'IESG=',IESG,'NU=',NUWRITE(2,*)'NODEXY‘DO10I=1,NWREAD(1,*)NODP_X(I),NODP_Y(I)WRITE(2,*)I,NODP_X(I),NODP_Y(I)CONTINUEWRITE(2,*)'ELEMENTHLHREF‘DO20NT=1,IESGREAD(1,*)BAR_HL(NT),BAR_HR(NT),BAR_EF(NT)WRITE(2,*)BAR_HL(NT),BAR_HR(NT),BAR_EF(NT)20CONTINUE!LISTDATA,NOWTOORGANIZETHEGLOBALMATRIX{輸入本體結構和性質的數(shù)據(jù)數(shù)據(jù)}4‘NW4‘IESG2‘NU01‘X;Y11100014100‘HI;HL;EF12100231001310010‘PX;PY10根據(jù)平面桁架的特性進一步展開為：{總體數(shù)據(jù)結構的設計}{輸入描述平面桁架本體結構和性質的數(shù)據(jù)}{輸出結構圖形}{計算各桿件的單元剛度矩陣}{形成結構總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項總外力向量}{回代求解，求出總位移向量}{打印位移}{計算各桿內力}{輸出各桿內力}{計算各桿件的單元剛度矩陣}相應位移向量SUBROUTINESTIF(NT,L,C,B,RD)USEPDATAINTEGERNT,I,JREALL,RD,C(4,4),B(4),CA,SA,DX,DYI=BAR_HL(NT)J=BAR_HR(NT) !LEFTANDRIGHTNODENUMBERDX=NODP_X(J)-NODP_X(I)DY=NODP_Y(J)-NODP_Y(I)L=SQRT(DX*DX+DY*DY)CA=DX/L;SA=DY/L !(*DIRECTINCOSINE*)B(1)=-CAB(2)=-SAB(3)=CAB(4)=SA!(*THEDISPLACEMENT->DEFORMATIONTRANSFORMATIONMATRIX*)RD=BAR_EF(NT)/L!(*BARSTIFFNESS*)DO10I=1,4DO10J=1,4C(I,J)=RD*B(I)*B(J)CONTINUE!(*CISTHESTIFFNESSMATRIXINGLOBALCOORDINATE*)RETURNEND !(*STIF*){計算各桿件的單元剛度矩陣}根據(jù)平面桁架的特性進一步展開為：{總體數(shù)據(jù)結構的設計}{輸入描述平面桁架本體結構和性質的數(shù)據(jù)}{輸出結構圖形}{計算各桿件的單元剛度矩陣}{形成結構總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項總外力向量}{回代求解，求出總位移向量}{打印位移}{計算各桿內力}{輸出各桿內力}{形成結構總剛度矩陣}{對總剛度矩陣R清零}{FORNT:=1~IESGDO

注：對所有桿件循環(huán)BEGIN調出STIF(NT,T,C,B,RD)；計算桿件在全局坐標系下單元剛度矩陣C}調出I0J0(NT,II)；計算對號入坐信息數(shù)組II}調出ASSEMB2(IB,C,R,NDISP,BV,V)；注：結構總剛度矩陣R已形成END;}

!確定這根桿件兩端點（HL端和HR端）位移在總位移向量中的編號SUBROUTINEI0J0(NT,II)USEPDATA

INTEGERNT,II(4),I,I0,J0I=BAR_HL(NT)-1I0=I+I!NOTETHENUMBERINGFORMULADO10I=1,210

II(I)=I0+I

I=BAR_HR(NT)-1J0=I+IDO20I=3,4II(I)=J0+I-2RETURNEND!(*I0J0*)

ui在總位移向量中是第2*（i-1）+1號uj在總位移向量中是第2*（j-1）+1號SUBROUTINESTIF(NT,L,C,B,RD)USEPDATAINTEGERNT,I,JREALL,RD,C(4,4),B(4),CA,SA,DX,DYI=BAR_HL(NT)J=BAR_HR(NT) !LEFTANDRIGHTNODENUMBERDX=NODP_X(J)-NODP_X(I)DY=NODP_Y(J)-NODP_Y(I)L=SQRT(DX*DX+DY*DY)CA=DX/L;SA=DY/L !(*DIRECTINCOSINE*)B(1)=-CAB(2)=-SAB(3)=CAB(4)=SA!(*THEDISPLACEMENT->DEFORMATIONTRANSFORMATIONMATRIX*)RD=BAR_EF(NT)/L!(*BARSTIFFNESS*)DO10I=1,4DO10J=1,4C(I,J)=RD*B(I)*B(J)CONTINUE!(*CISTHESTIFFNESSMATRIXINGLOBALCOORDINATE*)RETURNEND !(*STIF*)13581011182225330112106327對角元地址數(shù)組：半帶寬數(shù)組：3.2

-1.04.3

0.82.0

4.00.18.8

4.14.0

2.1

-0.81.20.30002.8

1.5003.1

0.703.7

-0.200-1.20-2.14.17.2剛度矩陣（二維存儲）：3.24.30.82.04.00.18.84.14.02.11.20.300002.81.5003.10.703.70004.17.2剛度矩陣（一維存儲）：存儲每一行的對角元素在一維存儲數(shù)組中R中的位置，存放在數(shù)組

V[BV:BV+NDISP]對角元素在一維數(shù)組R中的地址為V[BV+I]，而這個元素的值為

R[V[BV+I]]1358101119232634對角元地址數(shù)組：3.24.30.82.04.00.18.84.14.02.11.20.300002.81.5003.10.703.70004.17.2一維剛度數(shù)組：根據(jù)V可以求出：（1）第i行主元在R中的位置為V[BV+i]，數(shù)值為R[V[BV+i]]。如第3行，BV=0，i=3，V[3]=6，其值R[6]=r33。（2）第i行半帶寬（不計對角元）為：BDW[i]=V[BV+i]-V[BV+i-1]-1（3）在方陣中處于i行j列的元素rij在R中的地址為：ix=V[BV+i]-i+j當i>=j但是要求ix>V[BV+i-1]（i>=j），否則該元素為零，在R中并沒有存放。（4）第i行第一個非零元素所在列號為i1=V[BV+i-1]-V[BV+i]+i+1式中當i=1時V[BV]定義為零。（5）第i行對角元地址可用下面的遞推公式計算：V[i]=V[i-1]+BDW[i]+1式中BDW[i]為第i行的半帶寬。半帶寬選大的自然語言表示如下：{未知數(shù)總數(shù)NDISP已知}{安排數(shù)組V的存區(qū)}{半帶寬存區(qū)清零，準備選大}{對桿件循環(huán)}{調出I0J0（NT,II）}{為計算NT號桿對半帶寬作準備}{選出半帶寬}根據(jù)自然語言表示，半帶寬選大的程序段是:BV:=1;FORI:=0TONDISPDOV[BV+I]:=0注：半帶寬選大之前，對V數(shù)組清零FORNT:=1TOIESGDOBEGINI0J0(NT,II)桿位移對號FORI:=1TOIBDO!IB為數(shù)組C的下標界偶，即單元位移向量的元素數(shù)BEGINFORJ:=1TOIDO對桿件的行循環(huán)和列循環(huán)IA:=II[I];JA:=II[J];轉換到總剛度矩陣的行、列號IA與JAIFIA<JATHENBEGINK:=IA;IA:=JA;JA:=K;!保證總剛度矩陣元素的行號IA大于列號JAEND;IF(JA>0)AND(IA<=NDISP)THEN防止出界BEGINIFV[BV+IA]<IA-JATHENV[BV+IA]:=IA-JA;!選大END;END;END;{半帶寬已選出}在此基礎上就可以編制DIAGADR(NDISP,V)形成總剛度矩陣的對角元地址過程。功能：在給定總剛度矩陣R存放的基址IV，半帶寬數(shù)組V[BV:BV+NDISP]的條件下，求出總剛度矩陣各行的對角元地址，并仍就存放在V[BV+1]~V[BV+NDISP]中。 SUBROUTINEASSEMB2(IB,II,C,R,NDISP,V,PHASE) USEPDATA INTEGERIB,NDISP,PHASE,I,J,K,IA,JA,II(4),V(120) REALC(4,4),R(1200) DO100I=1,IB DO50J=1,I IA=II(I) JA=II(J) !TAKEOUTDISPNUMBER IF(IA.LT.JA)THEN K=IA IA=JA JA=K ENDIF IF((JA.GT.0).AND.(IA.LE.NDISP))THEN IF(PHASE.EQ.2)THEN K=V(IA)-IA+JA!(*THEADESSINR*) R(K)=R(K)+C(I,J)!(*ACCUMULATING*) ELSE IF(V(IA).LT.(IA-JA))V(IA)=IA-JA ENDIF ENDIF50CONTINUE100CONTINUE RETURN END!ASSEMB2 SUBROUTINEDIAGADR(NDISP,V)USEPDATAINTEGERV(120),IA!FROMSEMI-BANDWIDTHTODIAGONALADRESSOFGLOBALMATRIX!IV--BASEADRESSOFGLOBALMATRIX!V[BV..BV+NDISP]SEMI-BANDINBUTRETURNWITHDIAGONALADRESS!NDISP:GLOBALVARIABLEDO10IA=1,NDISP10V(IA)=V(IA)+V(IA-1)+1RETURNEND!DIAGADR{總剛度矩陣三角化}修改的平方根法運算工作量小且可以處理滿陣存儲的系數(shù)矩陣，可以使用較小的存儲空間。每增加一個右端項時只須增加一次回代求解，而無須重新分解系數(shù)矩陣。具體實施在形成結構總剛度矩陣之后。!NDISP—THEORDEROFMATRIXTOBELDLTED!V[BV..BV+NDISP]—THEARRAYOFDIAGONALADREDD,THEMATRIXTOBE!LDLTEDISINR[V[BV]+1..V[BV+NDISP]];!T[1..BDW]—THEWORKINGARRAY;BDWISTHEMAX-BANDOFRSUBROUTINELDLT1(NDISP,V,R,T)根據(jù)平面桁架的特性進一步展開為：{總體數(shù)據(jù)結構的設計}{輸入描述平面桁架本體結構和性質的數(shù)據(jù)}{輸出結構圖形}{計算各桿件的單元剛度矩陣}{形成結構總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項總外力向量}{回代求解，求出總位移向量}{打印位移}{計算各桿內力}{輸出各桿內力}{輸入荷載數(shù)據(jù)，形成右端項總外力向量}WRITE(2,*)"NODEDPXDPY"DO70I=1,NU J=I+I-2 DP(J+1)=0.0 READ(1,*)DP(J+1),DP(J+2) WRITE(2,*)DP(J+1),DP(J+2)CONTINUE!THELOADONNODESHAVEBEENREADIN4‘NW4‘IESG2‘NU01‘X;Y11100014100‘HI;HL;EF12100231001310010‘PX;PY10{回代求解，求出總位移向量}調用SLVEQ1(NDISP,V,R,DP)計算結構的總位移向量。

輸出節(jié)點位移：WRITE(2,*) 'NODEXY