平面桁架程序的計(jì)算原理及程序編制

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

角是桿件軸線與方向的夾角，由正向逆時(shí)針向轉(zhuǎn)至桿軸的角度為正。進(jìn)一步有內(nèi)力和節(jié)點(diǎn)位移之間的關(guān)系對(duì)于每根桿件（以兩端節(jié)點(diǎn)編號(hào)A和B定出角）應(yīng)用上述公式，有應(yīng)用上述各桿內(nèi)力和位移關(guān)系后，便可建立以位移為未知數(shù)的節(jié)點(diǎn)平衡方程。平面桁架結(jié)構(gòu)節(jié)點(diǎn)平衡方程為

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

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

寫成列向量{V}，系數(shù)排列成行向量{B}T，上式可以寫成如下形式其中把⊿L代入上面的公式，得到此時(shí)，已用全局坐標(biāo)系中的位移表達(dá)出桿的內(nèi)力N。由于最終得到的平衡方程都是相對(duì)于全局坐標(biāo)系建立的，上面所計(jì)算出的內(nèi)力與桿軸方向一致。假設(shè)在桿端作用平行于桿軸的力FA、FB，則由平衡方程FA和FB在x軸、y軸方向的分量分別為：把、、、排列成列向量，則有由此令得到全局坐標(biāo)系下的單元?jiǎng)偠染仃?-8結(jié)構(gòu)總剛度矩陣的形成假設(shè)平衡方程形式為展開上式可以寫成式中nn=2*NU；P1，P2，…，Pnn為加在節(jié)點(diǎn)上的外力。5-9桿件內(nèi)力的計(jì)算對(duì)任意一根桿件，如果桿的兩端位移已知。這些位移在總位移向量中i0+1，i0+2，j0+1，j0+2的位置上，表示如下進(jìn)而用下面公式計(jì)算結(jié)構(gòu)內(nèi)力平面桁架計(jì)算程序的編制矩陣位移法求解一般步驟：（1）結(jié)構(gòu)計(jì)算簡(jiǎn)圖。（節(jié)點(diǎn)、單元編號(hào)，建立一個(gè)統(tǒng)一的坐標(biāo)系等）（2）分析節(jié)點(diǎn)位移的力學(xué)特性，確定位移未知數(shù)。（3）建立每根桿件兩端位移和內(nèi)力的關(guān)系。（單根桿件的剛度矩陣）（4）根據(jù)每根桿件上的上述關(guān)系建立結(jié)構(gòu)可動(dòng)節(jié)點(diǎn)的平衡方程。（結(jié)構(gòu)總剛度矩陣）（5）求解平衡方程，得到節(jié)點(diǎn)位移。（6）根據(jù)求得的位移，利用每根桿件位移與內(nèi)力關(guān)系計(jì)算各桿的內(nèi)力。程序流程平面桁架計(jì)算程序的總功能表示{依據(jù)平面桁架及計(jì)算模型的有關(guān)數(shù)據(jù)，進(jìn)行結(jié)構(gòu)靜力計(jì)算，并輸出計(jì)算結(jié)果}根據(jù)平面桁架的特性進(jìn)一步展開為：{總體數(shù)據(jù)結(jié)構(gòu)的設(shè)計(jì)}{輸入描述平面桁架本體結(jié)構(gòu)和性質(zhì)的數(shù)據(jù)}{輸出結(jié)構(gòu)圖形}{計(jì)算各桿件的單元?jiǎng)偠染仃噠{形成結(jié)構(gòu)總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項(xiàng)總外力向量}{回代求解，求出總位移向量}{打印位移}{計(jì)算各桿內(nèi)力}{輸出各桿內(nèi)力}根據(jù)平面桁架的特性進(jìn)一步展開為：{總體數(shù)據(jù)結(jié)構(gòu)的設(shè)計(jì)}{輸入描述平面桁架本體結(jié)構(gòu)和性質(zhì)的數(shù)據(jù)}{輸出結(jié)構(gòu)圖形}{計(jì)算各桿件的單元?jiǎng)偠染仃噠{形成結(jié)構(gòu)總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項(xiàng)總外力向量}{回代求解，求出總位移向量}{打印位移}{計(jì)算各桿內(nèi)力}{輸出各桿內(nèi)力}{總體數(shù)據(jù)結(jié)構(gòu)的設(shè)計(jì)}平面桁架程序的總體數(shù)據(jù)結(jié)構(gòu)設(shè)計(jì)，包括全局量（標(biāo)識(shí)符）的存儲(chǔ)安排與說明。一個(gè)平面桁架結(jié)構(gòu)的數(shù)據(jù)有三個(gè)方面：結(jié)構(gòu)本體描述數(shù)據(jù)（NW,IESG,NU,X,Y,HL,HR）性質(zhì)數(shù)據(jù)（EF）荷載數(shù)據(jù)（PX,PY）NW為節(jié)點(diǎn)總數(shù)IESG桿件總數(shù)NU可動(dòng)節(jié)點(diǎn)總數(shù)X,Y節(jié)點(diǎn)坐標(biāo)HL,HR每根桿件兩端節(jié)點(diǎn)編號(hào)EF性質(zhì)數(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ù)結(jié)構(gòu)的設(shè)計(jì)}4‘NW4‘IESG2‘NU01‘X;Y11100014100‘HI;HL;EF12100231001310010‘PX;PY10節(jié)點(diǎn)總數(shù)：NW=4可動(dòng)節(jié)點(diǎn)數(shù)：NU=2位移未知數(shù)總數(shù)：NDISP=2*NU=4平面桁架結(jié)構(gòu)1342{輸入描述平面桁架本體結(jié)構(gòu)和性質(zhì)的數(shù)據(jù)}結(jié)構(gòu)本體描述數(shù)據(jù)（NW,IESG,NU,X,Y,HL,HR）性質(zhì)數(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{輸入本體結(jié)構(gòu)和性質(zhì)的數(shù)據(jù)數(shù)據(jù)}4‘NW4‘IESG2‘NU01‘X;Y11100014100‘HI;HL;EF12100231001310010‘PX;PY10根據(jù)平面桁架的特性進(jìn)一步展開為：{總體數(shù)據(jù)結(jié)構(gòu)的設(shè)計(jì)}{輸入描述平面桁架本體結(jié)構(gòu)和性質(zhì)的數(shù)據(jù)}{輸出結(jié)構(gòu)圖形}{計(jì)算各桿件的單元?jiǎng)偠染仃噠{形成結(jié)構(gòu)總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項(xiàng)總外力向量}{回代求解，求出總位移向量}{打印位移}{計(jì)算各桿內(nèi)力}{輸出各桿內(nèi)力}{計(jì)算各桿件的單元?jiǎng)偠染仃噠相應(yīng)位移向量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ì)算各桿件的單元?jiǎng)偠染仃噠根據(jù)平面桁架的特性進(jìn)一步展開為：{總體數(shù)據(jù)結(jié)構(gòu)的設(shè)計(jì)}{輸入描述平面桁架本體結(jié)構(gòu)和性質(zhì)的數(shù)據(jù)}{輸出結(jié)構(gòu)圖形}{計(jì)算各桿件的單元?jiǎng)偠染仃噠{形成結(jié)構(gòu)總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項(xiàng)總外力向量}{回代求解，求出總位移向量}{打印位移}{計(jì)算各桿內(nèi)力}{輸出各桿內(nèi)力}{形成結(jié)構(gòu)總剛度矩陣}{對(duì)總剛度矩陣R清零}{FORNT:=1~IESGDO

注：對(duì)所有桿件循環(huán)BEGIN調(diào)出STIF(NT,T,C,B,RD)；計(jì)算桿件在全局坐標(biāo)系下單元?jiǎng)偠染仃嘋}調(diào)出I0J0(NT,II)；計(jì)算對(duì)號(hào)入坐信息數(shù)組II}調(diào)出ASSEMB2(IB,C,R,NDISP,BV,V)；注：結(jié)構(gòu)總剛度矩陣R已形成END;}

!確定這根桿件兩端點(diǎn)（HL端和HR端）位移在總位移向量中的編號(hào)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號(hào)uj在總位移向量中是第2*（j-1）+1號(hào)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對(duì)角元地址數(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剛度矩陣（二維存儲(chǔ)）：3.24.30.82.04.00.18.84.14.02.11.20.300002.81.5003.10.703.70004.17.2剛度矩陣（一維存儲(chǔ)）：存儲(chǔ)每一行的對(duì)角元素在一維存儲(chǔ)數(shù)組中R中的位置，存放在數(shù)組

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

R[V[BV+I]]1358101119232634對(duì)角元地址數(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行半帶寬（不計(jì)對(duì)角元）為：BDW[i]=V[BV+i]-V[BV+i-1]-1（3）在方陣中處于i行j列的元素rij在R中的地址為：ix=V[BV+i]-i+j當(dāng)i>=j但是要求ix>V[BV+i-1]（i>=j），否則該元素為零，在R中并沒有存放。（4）第i行第一個(gè)非零元素所在列號(hào)為i1=V[BV+i-1]-V[BV+i]+i+1式中當(dāng)i=1時(shí)V[BV]定義為零。（5）第i行對(duì)角元地址可用下面的遞推公式計(jì)算：V[i]=V[i-1]+BDW[i]+1式中BDW[i]為第i行的半帶寬。半帶寬選大的自然語言表示如下：{未知數(shù)總數(shù)NDISP已知}{安排數(shù)組V的存區(qū)}{半帶寬存區(qū)清零，準(zhǔn)備選大}{對(duì)桿件循環(huán)}{調(diào)出I0J0（NT,II）}{為計(jì)算NT號(hào)桿對(duì)半帶寬作準(zhǔn)備}{選出半帶寬}根據(jù)自然語言表示，半帶寬選大的程序段是:BV:=1;FORI:=0TONDISPDOV[BV+I]:=0注：半帶寬選大之前，對(duì)V數(shù)組清零FORNT:=1TOIESGDOBEGINI0J0(NT,II)桿位移對(duì)號(hào)FORI:=1TOIBDO!IB為數(shù)組C的下標(biāo)界偶，即單元位移向量的元素?cái)?shù)BEGINFORJ:=1TOIDO對(duì)桿件的行循環(huán)和列循環(huán)IA:=II[I];JA:=II[J];轉(zhuǎn)換到總剛度矩陣的行、列號(hào)IA與JAIFIA<JATHENBEGINK:=IA;IA:=JA;JA:=K;!保證總剛度矩陣元素的行號(hào)IA大于列號(hào)JAEND;IF(JA>0)AND(IA<=NDISP)THEN防止出界BEGINIFV[BV+IA]<IA-JATHENV[BV+IA]:=IA-JA;!選大END;END;END;{半帶寬已選出}在此基礎(chǔ)上就可以編制DIAGADR(NDISP,V)形成總剛度矩陣的對(duì)角元地址過程。功能：在給定總剛度矩陣R存放的基址IV，半帶寬數(shù)組V[BV:BV+NDISP]的條件下，求出總剛度矩陣各行的對(duì)角元地址，并仍就存放在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{總剛度矩陣三角化}修改的平方根法運(yùn)算工作量小且可以處理滿陣存儲(chǔ)的系數(shù)矩陣，可以使用較小的存儲(chǔ)空間。每增加一個(gè)右端項(xiàng)時(shí)只須增加一次回代求解，而無須重新分解系數(shù)矩陣。具體實(shí)施在形成結(jié)構(gòu)總剛度矩陣之后。!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ù)平面桁架的特性進(jìn)一步展開為：{總體數(shù)據(jù)結(jié)構(gòu)的設(shè)計(jì)}{輸入描述平面桁架本體結(jié)構(gòu)和性質(zhì)的數(shù)據(jù)}{輸出結(jié)構(gòu)圖形}{計(jì)算各桿件的單元?jiǎng)偠染仃噠{形成結(jié)構(gòu)總剛度矩陣}{總剛度矩陣三角化}{輸入荷載數(shù)據(jù)，形成右端項(xiàng)總外力向量}{回代求解，求出總位移向量}{打印位移}{計(jì)算各桿內(nèi)力}{輸出各桿內(nèi)力}{輸入荷載數(shù)據(jù)，形成右端項(xiàng)總外力向量}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{回代求解，求出總位移向量}調(diào)用SLVEQ1(NDISP,V,R,DP)計(jì)算結(jié)構(gòu)的總位移向量。

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