C Ping Mian Gang Jia Cheng Xu-----GJ2001.FOR 2004.4.2 CHARACTER*12 FA,FB INTEGER MN(2,40),R(3,40),A(120),CM(6),U,V,W,S,T,T1,G,PI,PJ,PN REAL KM(6,6),C(6,6),CK(6,6),X(40),Y(40),EA(40),EI(40),FO(6), * K(4000),P1(4,40),P2(4,100),P3(7,40),PP(120),CN(6),CV(6),AL COMMON MI,MJ,PI,PJ,PN C ㈠ 分别输入数据文件名,计算结果文件名 WRITE(*,*) 'Please input data file name:' c read datafile--0020 READ(*,'(A)') FA OPEN (20,FILE=FA,STATUS='OLD') WRITE(*,*) 'Please input result file name:' c read outfile--0021 READ(*,'(A)') FB OPEN (21,FILE=FB,STATUS='NEW') C ㈡ 输入原始数据，形成结点位移编号 WRITE (21,20) 20 FORMAT (1X,'Program name: GJ2001.FOR---2004.4.2') WRITE (21,30) FA 30 FORMAT (1X,'Data file name: ',A) WRITE (21,40) FB 40 FORMAT (1X,'Result file name: ',A) c read data--1 READ (20,*) MI,MJ,PI,PJ,PN WRITE(21,50) MI,MJ,PI,PJ,PN 50 FORMAT(1X,'MI=',I3,' MJ=',I3,' PI=',I3,' PJ=',I3,' PN=',I3) c read data--2 READ (20,*) ((R(I,J),I=1,3),J=1,MJ) WRITE(21,60) ((R(I,J),I=1,3),J=1,MJ) 60 FORMAT(1X,'R(i,j)='/8X,50(/1X,12I6)) c read data--3 READ (20,*) ((MN(I,J),I=1,2),J=1,MI) WRITE(21,70) ((MN(I,J),I=1,2),J=1,MI) 70 FORMAT(1X,'MN(i,j)='/8X,50(/1X,10I6)) c read data--4 READ (20,*) (X(I),Y(I),I=1,MJ) WRITE(21,80) (X(I),Y(I),I=1,MJ) 80 FORMAT(1X,'X(i)--Y(i)='/1X,50(/1X,10F7.3)) c read data--5 READ (20,*) (EA(I),EI(I),I=1,MI) WRITE(21,90) (EA(I),EI(I),I=1,MI) 90 FORMAT(1X,'EA(i)--EI(i)='/1X,50(/1X,6E12.3)) IF (PJ.NE.0) THEN c read data--6 READ (20,*)((P1(J,I),J=1,4),I=1,PJ) WRITE(21,100) ((P1(J,I),J=1,4),I=1,PJ) 100 FORMAT(1X,'P1(i,j)='/1X,50(/1X,F4.1,3E13.3)) END IF IF (PN.NE.0) THEN c read data--7 READ (20,*) ((P2(J,I),J=1,4),I=1,PN) WRITE(21,110) ((P2(J,I),J=1,4),I=1,PN) 110 FORMAT(1X,'P2(i,j)='/1X,50(/1X,F4.1,3E13.3)) END IF N=0 MAX=0 DO 150 I=1,MJ DO 150 J=1,3 J1=R(J,I) IF (J1.EQ.0) GOTO 150 IF (J1.EQ.1) THEN N=N+1 R(J,I)=N ELSE IF (J1.GT.1) THEN S=J1 T=R(J,S) R(J,I)=T ELSE END IF WRITE (*,140) I,J,R(J,I) 140 FORMAT (1X,'R(',I3,',',I3,')=',I3) 150 CONTINUE WRITE (*,160)N 160 FORMAT (1X,'N=',I3) DO 170 I=1,N A(I)=0 170 CONTINUE DO 200 M=1,MI CALL DJWB(M,MN,R,CM) DO 190 S=1,6 I1=CM(S) IF (I1.EQ.0) GOTO 190 DO 180 T=1,6 W=CM(T) IF (W.EQ.0.OR.W.GT.I1) GOTO 180 U=I1-W+1 IF (U.GT.A(I1)) A(I1)=U 180 CONTINUE 190 CONTINUE 200 CONTINUE A(1)=1 DO 210 I=2,N IF (A(I).GT.MAX) MAX=A(I) A(I)=A(I)+A(I-1) 210 CONTINUE WRITE (*,220)(I,A(I),I=1,N) 220 FORMAT(1X,4('A(',I3,')=',I3,2X)) WRITE (*,230) MAX 230 FORMAT (4X,'MAX=',I3) C ㈢ 形成结构的刚度矩阵［Ｋ］ N1=A(N) DO 240 I=1,N1 K(I)=0.0 240 CONTINUE DO 290 M=1,MI CALL DCS(M,MN,X,Y,AL,CO,SI) CALL DGJ1 (M,AL,EA,EI,KM) IF (SI.EQ.0) GOTO 245 CALL BH (CO,SI,C) CALL DGJ2 (C,CK,KM) 245 CALL DJWB (M,MN,R,CM) WRITE (*,250)M 250 FORMAT (4X,'KM(I,J)=(',I3,')') WRITE (*,260)((KM(S,T),S=1,6),T=1,6) 260 FORMAT (1X,6E12.4) DO 280 S=1,6 IF (CM(S).EQ.0) GOTO 280 W=CM(S) DO 270 T=1,S IF (CM(T).EQ.0) GOTO 270 G=CM(T) IF (W.GE.G) U=A(W)-W+G IF (W.LT.G) U=A(G)-G+W K(U)=K(U)+KM(S,T) 270 CONTINUE 280 CONTINUE 290 CONTINUE WRITE (*,*)'K(N)=' WRITE (*,300)(K(S),S=1,N1) 300 FORMAT (1X,5E12.4) C ㈣ 形成结构的结点荷载向量｛Ｐ｝ DO 310 I=1,N PP(I)=0.0 310 CONTINUE IF (PJ.EQ.0) GOTO 350 WRITE (*,*)'P1(I,J)=' WRITE (*,320)((P1(I,J),I=1,4),J=1,PJ) 320 FORMAT (1X,F4.1,3F10.4) DO 340 I=1,PJ S=IFIX(P1(1,I)) DO 330 J=1,3 T=J+1 U=R(J,S) IF(U.EQ.0) GOTO 330 PP(U)=PP(U)+P1(T,I) 330 CONTINUE 340 CONTINUE 350 IF(PN.EQ.0) GOTO 410 DO 370 I=1, PI DO 360 J=1,6 P3(J,I)=0.0 360 CONTINUE 370 CONTINUE CALL DGL(P2,P3,X,Y,MN,EA,EI,KM,FO,CN,I1,J1,AL,PN) DO 400 I=1,PI W=IFIX(P3(7,I)) CALL DCS (W,MN,X,Y,AL,CO,SI) CALL BH (CO,SI,C) CALL DJWB (W,MN,R,CM) DO 390 S=1,6 CV(S)=0.0 U=CM(S) IF (U.EQ.0) GOTO 390 DO 380 T=1,6 CV(S)=CV(S)+C(T,S)*P3(T,I) 380 CONTINUE PP(U)=PP(U)-CV(S) 390 CONTINUE 400 CONTINUE C ㈤ 三角分解法解方程，求结构的结点位移向量｛Δ｝ T=0 410 DO 440 I=1,N S=A(I) IF(I.EQ.1) GOTO 412 T=A(I-1) 412 V=I-S+T+1 DO 430 J=V,I J2=A(J) T1=A(J-1) W=S-I+J J1=J-1 Z=0.0 DO 420 U=V,J1 I1=J-J2+T1 IF (I1.GE.U) GOTO 420 I2=A(U) G=S-I+U IJ=J2-J+U Z=Z+K(G)*K(IJ)/K(I2) 420 CONTINUE K(W)=K(W)-Z 430 CONTINUE 440 CONTINUE T=0 DO 460 I=1,N S=A(I) IF(I.EQ.1) GOTO 442 T=A(I-1) 442 V=I-S+T+1 I1=I-1 Z=0.0 DO 450 J=V,I1 W=S-I+J Z=Z+K(W)*PP(J) 450 CONTINUE PP(I)=(PP(I)-Z)/K(S) 460 CONTINUE DO 480 I=N,1,-1 U=A(I) I1=I+1 Z=0.0 DO 470 J=I1,N S=A(J) T=A(J-1) IH=S-T W=J-I+1 G=S-J+I IF (W.GT.IH) GOTO 470 Z=Z+K(G)*PP(J)/K(U) 470 CONTINUE PP(I)=PP(I)-Z 480 CONTINUE WRITE (*,490) WRITE (21,490) 490 FORMAT (7X,'I',7X,'NX',13X,'NY',13X,'NF'/) DO 520 I=1,MJ DO 500 J=1,3 CN(J)=0.0 T=R(J,I) IF (T.EQ.0) GOTO 500 CN(J)=PP(T) 500 CONTINUE WRITE (*,510) I,CN(1),CN(2),CN(3) WRITE (21,510)I,CN(1),CN(2),CN(3) 510 FORMAT (2X,I6,3E13.4) 520 CONTINUE C ㈥ 求单元的杆端内力并输出｛｝ WRITE (*,530) WRITE (21,530) 530 FORMAT (3X,'I',6X,'NI',10X,'QI',10X,'MI', * 11X,'NJ',10X,'QJ',10X,'MJ'/) DO 640 I=1,MI CALL DCS (I,MN,X,Y,AL,CO,SI) CALL DGJ1 (I,AL,EA,EI,KM) CALL DJWB (I,MN,R,CM) DO 540 J=1,6 IF (CM(J).NE.0) THEN U=CM(J) D=PP(U) CN(J)=D ELSE CN(J)=0.0 END IF 540 CONTINUE CALL BH (CO,SI,C) DO 560 J=1,6 CV(J)=0.0 DO 550 S=1,6 CV(J)=CV(J)+C(J,S)*CN(S) 550 CONTINUE 560 CONTINUE DO 580 J=1,6 FO(J)=0.0 DO 570 S=1,6 FO(J)=FO(J)+KM(J,S)*CV(S) 570 CONTINUE 580 CONTINUE DO 590 J=1,PI M=P3(7,J) IF (M.EQ.I) GOTO 600 590 CONTINUE GOTO 620 600 DO 610 S=1,6 FO(S)=FO(S)+P3(S,J) 610 CONTINUE 620 FO(1)=-FO(1) FO(3)=-FO(3) FO(5)=-FO(5) FO(6)=-FO(6) WRITE (*,630) I,(FO(J),J=1,6) WRITE (21,630)I,(FO(J),J=1,6) 630 FORMAT (1X,I2,6F12.3) 640 CONTINUE CLOSE (20) END C ㈦ 子程序⑴形成杆件坐标的单元刚度矩阵［］ SUBROUTINE DGJ1(M,AL,EA,EI,KM) COMMON MI,MJ DIMENSION EA(MI),EI(MI) INTEGER S,T REAL AL,B1,B2,B3,B4,KM(6,6) B1=EA(M)/AL B2=EI(M)/AL B3=6.0*B2/AL B4=2.0*B3/AL DO 10 S=1,6 DO 20 T=1,6 KM(S,T)=0.0 20 CONTINUE 10 CONTINUE KM(1,1)=B1 KM(1,4)=-B1 KM(2,2)=B4 KM(2,3)=B3 KM(2,5)=-B4 KM(2,6)=B3 KM(3,3)=4.0*B2 KM(3,5)=-B3 KM(3,6)=2.0*B2 KM(4,4)=B1 KM(5,5)=B4 KM(5,6)=-B3 KM(6,6)=4.0*B2 DO 30 S=2,6 DO 40 T=1,S-1 KM(S,T)=KM(T,S) 40 CONTINUE 30 CONTINUE END C ㈧ 子程序⑵形成坐标变换矩阵［Ｃ］ SUBROUTINE BH (CO,SI,C) INTEGER S,T REAL CO,SI,C(6,6) DO 10 S=1,6 DO 20 T=1,6 C(S,T)=0.0 20 CONTINUE 10 CONTINUE C(1,1)=CO C(1,2)=SI C(2,1)=-SI C(2,2)=CO C(3,3)=1.0 DO 30 S=1,3 DO 40 T=1,3 C(S+3,T+3)=C(S,T) 40 CONTINUE 30 CONTINUE END C ㈨ 子程序⑶形成结构坐标的单元刚度矩阵［ｋ］ SUBROUTINE DGJ2(C,CK,KM) DIMENSION C(6,6),CK(6,6) INTEGER S,T REAL KM(6,6) DO 10 S=1,6 DO 20 T=1,6 CK(S,T)=0.0 DO 30 M=1,6 CK(S,T)=CK(S,T)+C(M,S)*KM(M,T) 30 CONTINUE 20 CONTINUE 10 CONTINUE DO 40 S=1,6 DO 50 T=1,6 KM(S,T)=0.0 DO 60 M=1,6 KM(S,T)=KM(S,T)+CK(S,M)*C(M,T) 60 CONTINUE 50 CONTINUE 40 CONTINUE END C ㈩ 子程序⑷形成单元定位向量 SUBROUTINE DJWB(I,MN,R,CM) COMMON MI,MJ DIMENSION MN(2,MI),R(3,MJ),CM(6) INTEGER R,CM I1=MN(1,I) I2=MN(2,I) DO 10 J=1,3 J1=R(J,I1) CM(J)=J1 J1=R(J,I2) CM(J+3)=J1 10 CONTINUE END C （十一） 子程序⑸形成单元固端力向量｛。｝ SUBROUTINE DGL (P2,P3,X,Y,MN,EA,EI,KM,FO,CN,I1,J1,AL,PN) COMMON MI,MJ DIMENSION X(MJ),Y(MJ),MN(2,MI),EA(MI),EI(MI),P2(4,100),P3(7,50), * KM(6,6),FO(6),CN(6) INTEGER I1,J1,IL,MN,PN REAL KM,AL,L1,L2,XL,YL,Q,S,Q1,Q2,S1 I1=0 J1=0 DO 200 I2=1,PN M=P2(1,I2) IL=P2(2,I2) Q=P2(3,I2) S=P2(4,I2) IU=MN(1,M) IV=MN(2,M) XL=X(IV)-X(IU) YL=Y(IV)-Y(IU) AL=SQRT(XL*XL+YL*YL) Q1=Q/AL Q2=Q1*Q1 L1=AL-Q L2=L1/AL IF (M.EQ.J1) GOTO 20 I1=I1+1 J1=M P3(7,I1)=M 20 DO 30 J2=1,6 FO(J2)=0.0 30 CONTINUE GOTO (40,50,60,70,80,90,120,160) IL 40 FO(2)=-S*L2*L2*(1.0+2.0*Q1) FO(5)=-S*Q2*(1.0+2.0*L2) FO(3)=-S*Q*L2*L2 FO(6)=S*Q2*L1 GOTO 180 50 FO(2)=6.0*S*Q1*L2/AL FO(5)=-FO(2) FO(3)=S*L2*(3.0*Q1-1.0) FO(6)=S*Q1*(3.0*L2-1.0) GOTO 180 60 S1=0.5*S*Q FO(2)=-S1*(2.0-2.0*Q2+Q1*Q2) FO(5)=-S1*Q2*(2.0-Q1) S1=S1*Q/6.0 FO(6)=S1*(4.0*Q1-3.0*Q2) FO(3)=-S1*(6.0-8.0*Q1+3.0*Q2) GOTO 180 70 Q2=0.25*Q*S L1=Q1*Q1 L2=Q*Q1 FO(2)=-Q2*(2.0-3.0*L1+1.6*L1*Q1) FO(5)=-Q2*L1*(3.0-1.6*Q1) S1=Q*Q2/1.5 FO(3)=-S1*(2.0-3.0*Q1+1.2*L1) FO(6)=Q2*L2*(1.0-0.8*Q1) GOTO 180 80 FO(1)=-S*L2 FO(4)=-S*Q1 GOTO 180 90 L2=0.5*Q1 S1=1.0-L2 FO(1)=-S*S1*Q FO(4)=-0.5*S*Q*Q1 GOTO 180 c read data--8 120 READ (20,*)(CN(J),J=1,6) WRITE (21,130) M 130 FORMAT (1X,'CN(',I3,')=') WRITE (21,140) (CN(J),J=1,6) 140 FORMAT (1X,20(/1X,6E12.3)) IJ=M CALL DGJ1 (M,AL,EA,EI,KM) DO 150 IS=1,6 DO 150 JS=1,6 FO(IS)=FO(IS)+KM(IS,JS)*CN(JS) 150 CONTINUE GOTO 180 c read data--9 160 READ (20,*) T1,T2,T0 WRITE (21,170) T1,T2,T0 170 FORMAT (1X,'T1=',F7.2,' T2=',F7.2,' T0=',F7.2) FO(1)=EA(M)*Q*T0 FO(4)=-FO(1) FO(3)=Q*(T1-T2)*EI(M)/S FO(6)=-FO(3) 180 DO 190 J=1,6 P3(J,I1)=P3(J,I1)+FO(J) 190 CONTINUE 200 CONTINUE END C （十 二）子程序⑹形成单元常数 SUBROUTINE DCS (I,MN,X,Y,AL,CO,SI) COMMON MI,MJ DIMENSION MN(2,MI),X(MJ),Y(MJ) REAL AL,LX,LY,CO,SI I1=MN(1,I) I2=MN(2,I) LX=X(I2)-X(I1) LY=Y(I2)-Y(I1) AL=SQRT(LX*LX+LY*LY) CO=LX/AL SI=LY/AL END