【正文】 ,3),D(3) DOUBLE PRECISION A,D,B DO 20 J=2,N !對列循環(huán) JJ=ID(J) J1=J1 MJ=JJJ+ID(J1)+1 !第J列第一個非零元素的行號 IF() GOTO 20 !若MJJ1，則z的第J個值即為B的第J個值 DO 10 I=MJ,J1 JI=JJJ+I10 B(J)=B(J)A(JI)*B(I) !求z的第J個值，存放在B(N)數(shù)組。這里B(N)數(shù)組對應(yīng)主程序中的P(N)數(shù)組20 CONTINUE DO 30 I=1,N II=ID(I) B(I)=B(I)/A(II) !求y的第J個值，存放在B(N)數(shù)組30 CONTINUE NW=1 !半帶寬置初值為1 DO 35 I=2,N IW=ID(I)ID(I1) IF() NW=IW !計算最大列半帶寬35 CONTINUE N1=N1 DO 50 I=N1,1,1 !對行循環(huán) MIN=I+NW1 IF() MIN=N DO 40 J=I+1,MIN !對列循環(huán) JJ=ID(J) MJ=JJJ+ID(J1)+1 !第J列第一個非零元素的行號 IF() GOTO 40 !判斷是否為帶內(nèi)元素 JI=JJJ+I B(I)=B(I)A(JI)*B(J) !求x的第I個值，存放在B(N)數(shù)組。這里x的值即為結(jié)點位移40 CONTINUE50 CONTINUE WRITE(2,60)60 FORMAT(//2X,39。JOINT DISPLACEMENTS39。/5X,39。JOINT39。,12X,39。U39。,14X,39。V39。,11X,39。ceta39。) DO 80 I=1,NJ DO 70 J=1,3 D(J)= L=JN(I,J) !取位移號 IF() D(J)=B(L) !給位移賦值70 CONTINUE WRITE(2,75) I,D(1),D(2),D(3) !輸出第I號結(jié)點的位移值75 FORMAT(2X,I6,4X,)80 CONTINUE END SUBROUTINE MQN(M,NE,NJ,NAI,N,NF,E,X,Y,JE,MT,AI,JN,PF,P,F,H)!計算單元桿端力，存放在數(shù)組F(6)中 DIMENSION X(NJ),Y(NJ),JE(NE,2),AI(NAI,2),P(N),JN(NJ,3),PF(NF,4),JC(6),FO(6),D(6),FD(6),amp。F(6),KE(6,6),MT(NE),H(NE) DOUBLE PRECISION P,FO,D,FD,F,KE,BL,SI,CO CALL SCL(M,NE,NJ,X,Y,JE,BL,SI,CO) !求單元常數(shù) CALL ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE) !形成global坐標(biāo)下的單元剛度矩陣KE(6,6) CALL EJC(M,NE,NJ,JE,JN,JC) !形成單元定位向量JC DO 10 I=1,6 L=JC(I) !取位移編號 D(I)= IF() D(I)=P(L) !取global坐標(biāo)下的桿端位移10 CONTINUE DO 20 I=1,6 FD(I)= DO 30 J=1,6 FD(I)=FD(I)+KE(I,J)*D(J) !計算桿端力30 CONTINUE20 CONTINUE F(1)=FD(1)*CO+FD(2)*SI !將桿端力轉(zhuǎn)化到local坐標(biāo)下 F(2)=FD(1)*SI+FD(2)*CO !將桿端力轉(zhuǎn)化到local坐標(biāo)下 F(3)=FD(3) !將桿端力轉(zhuǎn)化到local坐標(biāo)下 F(4)=FD(4)*CO+FD(5)*SI !將桿端力轉(zhuǎn)化到local坐標(biāo)下 F(5)=FD(4)*SI+FD(5)*CO !將桿端力轉(zhuǎn)化到local坐標(biāo)下 F(6)=FD(6) !將桿端力轉(zhuǎn)化到local坐標(biāo)下 IF()THEN DO 40 I=1,NF L=INT(PF(I,1)) IF()THEN CALL EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H) DO 50 J=1,6 F(J)=F(J)+FO(J) !如果有非結(jié)點荷載，桿端力中要加入等效結(jié)點荷載50 CONTINUE ENDIF40 CONTINUE ENDIF END SUBROUTINE EJC(M,NE,NJ,JE,JN,JC) !形成單元定位向量JC，JC存放的是單元兩端點的位移號 DIMENSION JE(NE,2),JN(NJ,3),JC(6) I=JE(M,1) !取M單元起始結(jié)點號 J=JE(M,2) !取M單元終止結(jié)點號 DO 10 K=1,3 JC(K)=JN(I,K) !取單元位移號賦值給定位向量JC JC(K+3)=JN(J,K) !取單元位移號賦值給定位向量JC10 CONTINUE END SUBROUTINE ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE) !形成global坐標(biāo)下的單元剛度矩陣KE(6,6) DIMENSION MT(NE),AI(NAI,2),KE(6,6) DOUBLE PRECISION BL,SI,CO,KE I=MT(M) !取M單元截面類型號，MT——單元截面類型號 C1=E*AI(I,1)/BL C2=*E*AI(I,2)/BL C3=*C2/BL C4=*C3/BL S1=C1*CO*CO+C4*SI*SI S2=(C1C4)*SI*CO S3=C3*SI S4=C1*SI*SI+C4*CO*CO S5=C3*CO S6=C2 KE(1,1)=S1 KE(1,2)=S2 KE(1,3)=S3 KE(1,4)=S1 KE(1,5)=S2 KE(1,6)=S3 KE(2,2)=S4 KE(2,3)=S5 KE(2,4)=S2 KE(2,5)=S4 KE(2,6)=S5 KE(3,3)=*S6 KE(3,4)=S3 KE(3,5)=S5 KE(3,6)=S6 KE(4,4)=S1 KE(4,5)=S2 KE(4,6)=S3 KE(5,5)=S4 KE(5,6)=S5 KE(6,6)=*S6 DO 10 I=1,5 DO 10 J=I+1,610 KE(J,I)=KE(I,J) ! WRITE(2,20)((KE(I,J),J=1,6),I=1,6)!20 FORMAT(39。KE39。,3X,) END SUBROUTINE EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H)!將非結(jié)點荷載轉(zhuǎn)化為等效結(jié)點荷載存放在FO(6)中 DIMENSION PF(NF,4),MT(NE),AI(NAI,2),FO(6),H(NE) DOUBLE PRECISION BL,FO IF((PF(I,1)M).)GOTO 110 !判斷此非結(jié)點荷載是否為溫度荷載 NO=INT(PF(I,2)) !取此非結(jié)點荷載的類型號 Q=PF(I,3) C=PF(I,4) IF() GOTO 3 !判斷此非結(jié)點荷載是否為1－6號荷載 B=BLC C1=C/BL C2=C1*C1 C3=C1*C23 DO 5 J=1,65 FO(J)= GOTO(10,20,30,40,50,60,70,80,90,100),NO !對不同的非結(jié)點荷載分別計算其等效結(jié)點荷載10 FO(2)=Q*C*(+C3/) FO(3)=Q*C*C*(*C1/+*C2) FO(5)=Q*C*C2*(*C1) FO(6)=Q*C*C*C1*(*C1) RETURN20 FO(2)=Q*B*B*(+*C1)/BL/BL FO(3)=Q*C*B*B/BL/BL FO(5)=Q*C2*(+*B/BL) FO(6)=Q*C2*B RETURN 30 FO(2)=*Q*C1*B/BL/BL FO(3)=Q*B*(*B/BL)/BL FO(5)=*Q*C1*B/BL/BL FO(6)=Q*C1*(*C1) RETURN40 FO(2)=Q*C*(*C2+*C3) FO(3)=Q*C*C*(*C1+*C2) FO(5)=Q*C*C2*(*C1) FO(6)=Q*C*C*C1*(*C1) RETURN50 FO(1)=Q*C*(*C1) FO(4)=*Q*C*C1 RETURN60 FO(1)=Q*B/BL FO(4)=Q*C170 L=INT(C) K=MT(M) S=E*AI(K,1)*Q/BL FO(L)=S IF()FO(4)=S IF()FO(1)=S RETURN80 L=INT(C) K=MT(M) FO(L)=*E*AI(K,2)*Q/BL/BL/BL IF()FO(5)=FO(2) IF()FO(2)=FO(5) FO(3)=*BL*FO(2) FO(6)=FO(3) RETURN90 L=INT(C) K=MT(M) S=*E*AI(K,2)*Q/BL FO(L)=*S IF()FO(6)=S IF()FO(3)=S FO(2)=*S/BL FO(5)=FO(2) RETURN100 L=INT(C) K=MT(M) S=E*AI(K,1)*Q/BL FO(L)=S IF()FO(4)=S IF()FO(1)=S RETURN110 K=MT(M) FO(1)=*E*AI(K,1)*PF(I,4)*(PF(I,2)+PF(I,3)) FO(2)= FO(3)=E*AI(K,2)*PF(I,4)*(PF(I,2)+PF(I,3))/H(M) FO(4)=FO(1) FO(5)= FO(6)=FO(3) END**************************************************************************************算例： 如圖所示平面剛架，尺寸見圖。其中梁的截面積,。梁的截面積。支座沉降當(dāng)做單獨(dú)的工況計算。輸入數(shù)據(jù)文件（）：COMPOSITED PLANE FRAME STRUCTURE ANALYSISNE NJ NS NAI NL E6 7 3 2 2 ********************************************************COORDINATES OF JOINTJOINT X(I) Y(I)1 2 3 4 5 6 7 *********************************************************INFORMATION OF ELEMENTSELEMENT JOINTI JOINTJ AITYPE H1 2 4 1 02 1 3 2 03 4 5 2 04 3 5 1 05 3 6 2 06 5 7 2 0*******************************