平面刚架(Fortran77)电算程序

C 主程序

C (一)输入原始数据

DIMENSION JE(2,100),JN(3,100),JC(6),EA(100),EI(100),X(100),
* Y(100),PJ(2,50),PF(4,100)
REAL KE(6,6),KD(6,6),T(6,6),P(300),KB(300,20),F(6),F0(6),
* D(6),BL,SI,C0,S,C
OPEN(5,FILE='RPF.DAT')
READ(5,*)NE,NJ,N,NW,NPJ,NPF
READ(5,*)(X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ) !结点
READ(5,*)((JE(I,J),I=1,2),EA(J),EI(J),J=1,NE) !单元
IF(NPJ .NE. 0) READ(5,*) ((PJ(I,J),I=1,2),J=1,NPJ) !结点荷载
IF(NPF .NE. 0) READ(5,*) ((PF(I,J),I=1,4),J=1,NPF) !非结点荷载
WRITE(*,10) NE,NJ,N,NW,NPJ,NPF
WRITE(*,20) (J,X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ)
WRITE(*,30) (J,(JE(I,J),I=1,2),EA(J),EI(J),J=1,NE)
IF(NPJ .NE. 0) WRITE(*,40) ((PJ(I,J),I=1,2),J=1,NPJ)
IF(NPF .NE. O) WRITE(*,50) ((PF(I,J),I=1,4),J=1,NPF)
10 FORMAT(/6X,'NE=',I5,2X,'NJ=',I5,2X,'N=',I5,2X,'NW=',I5,2X,'NPJ=',
* I5,2X,'NPF=',I5)
20 FORMAT(/7X,'NODE',7X,'X',11X,'Y',12X,'XX',8X,'YY',8X,
* 'ZZ' /(1X,I10,2F12.4,3I10)) !结点
30 FORMAT(/4X,'ELEMENT',4X,'NODE-I',4X,'NODE-J',11X,
* 'EA',13X,'EI'/(1X,3I10,2E15.6)) !单元
40 FORMAT(/7X,'CODE',7X,'PX-PY-PM' /(1X,F10.0,F15.4))
50 FORMAT(/4X,'ELEMENT',7X,'IND',10X,'A',14X,'Q'/
* (1X,2F10.0,2F15.4))


C （二）形成总结点荷载向量
DO 55 I=1,N
55 P(I)=0.0
IF(NPJ .EQ. 0) GOTO 65
DO 60 I=I,NPJ
L=PJ(1,I)
60 P(L)=PJ(2,I) !结点荷载
65 IF(NPF .EQ. 0) GOTO 90
DO 70 I=1,NPF
M=PF(1,I) !M为单元编号
CALL SCL(M,NE,NJ,BL,SI,C0,JE,X,Y)
CALL EFX(I,NPF,BL,PF,F0)
CALL CTM(SI,C0,T)
CALL EJC(M,NE,NJ,JE,JN,JC)
DO 75 L=1,6
S=0.0
DO 80 K=1,6
80 S=S-T(K,L)*F0(K)
F(L)=S
75 CONTINUE !整体坐标系下的等效结点荷载
DO 85 J=1,6
L=JC(J)
IF(L .EQ. 0) GOTO 85
P(L)=P(L)+F(J)
85 CONTINUE
70 CONTINUE


C （三）形成整体刚度矩阵
90 DO 95 I=1,N
DO 100 J=1,NW
100 KB(I,J)=0.0
95 CONTINUE
DO 105 M=1,NE
CALL SCL(M,NE,NJ,BL,SI,C0,JE,X,Y)
CALL CTM(ST,C0,T)
CALL ESM(M,NE,BL,EA,EI,KD) !局部坐标下的单刚
CALL EJC(M,NE,NJ,JE,JN,JC)
DO 110 I=1,6
DO 115 J=1,6
S=0.0
DO 120 L=1,6
DO 125 K=1,6
125 S=S+T(L,I)*KD(L,K)*T(K,J)
120 CONTINUE
KE(I,J)=S
115 CONTINUE
110 CONTINUE !整体坐标下的单刚
DO 130 L=1,6
I=JC(L)
IF(I .EQ. 0) GOTO 130
DO 135 K=1,6
J=JC(K)
IF(J .EQ. 0 .OR. J .LT. I) GOTO 135
JJ=J-I+1
KB(I,JJ)=KB(I,JJ)+KE(L,K) !此时为半带宽存储，所以有此限制。
135 CONTINUE
130 CONTINUE
105 CONTINUE

C （四）解线性方程组
N1=N-1
DO 140 K=1,N1
IM=K+NW-1
IF(N .LT. IM) IM=N
I1=K+1
DO 145 I=I1,IM
L=I-K+1
C=KB(K,L)/KB(K,1)
JM=NW-L+1
DO 150 J=1,JM
JJ=J+I-K
150 KB(I,J)=KB(I,J)-C*KB(K,JJ)
145 P(I)=P(I)-C*P(K)
140 CONTINUE
P(N)=P(N)/KB(N,1)
DO 155 K=1,N1
I=N-K
JM=K+1
IF(NW .LT. JM) JM=NW
DO 160 J=2,JM
L=J+I-1
160 P(I)=P(I)-KB(I,J)*P(L)
155 P(I)=P(I)/KB(I,1)
WRITE(*,165)
165 FORMAT(/7X,'NODE',10X,'U',14X,'V',11X,'CETA')
DO 170 I=1,NJ
DO 175 J=1,3
D(J)=0.0
L=JN(J,I)
IF(L .EQ. 0) GOTO 175
D(J)=P(L)
175 CONTINUE
WRITE(*,180) I,D(1),D(2),D(3)
180 FORMAT(1X,I10,3E15.6)
170 CONTINUE

C （五）求单元杆端内力
WRITE(*,200)
200 FORMAT(/4X,' ELEMENT',13X, 'N', 17X, 'Q',17X,'M')
DO 205 M=1,NE
CALL SCL(M,NE,NJ,BL,SI,C0,JE,X,Y)
CALL ESM(M,NE,BL,EA,EI,KD)
CALL CTM(SI,C0,T)
DO 210 I=1,6
L=JC(I)
D(I)=0.0
IF(L .EQ. 0) GOTO 210
D(I)=P(L)
210 CONTINUE
DO 220 I=1,6
F(I)=0.0
DO 230 J=1,6
DO 240 K=1,6
240 F(I)=F(I)+KD(I,J)*T(J,K)*D(K)
230 CONTINUE
220 CONTINUE
IF(NPF .EQ. 0) GOTO 270
DO 250 I=1,NPF
L=PF(1,I)
IF(M .NE. L) GOTO 250
CALL EFX(I,NPF,BL,PF,F0)
DO 260 J=1,6
260 F(J)=F(J)+F0(J)
250 CONTINUE
270 WRITE(*,280) M,(F(I),I=1,6)
280 FORMAT(/1X,I10,3X,'N1=',F12.4,3X,'Q1=',F12.4,3X,'M1=',
* F12.4/14X,'N2=',F12.4,3X,'Q2=',F12.4,3X,'M2=',F12.4)
205 CONTINUE
CLOSE(5)
STOP
END

C （六）形成单元定位向量
SUBROUTINE EJC(M,NE,NJ,JE,JN,JC)
DIMENSION JE(2,NE),JN(3,NJ),JC(6)
J1=JE(1,M)
J2=JE(2,M)
DO 10 I=1,3
JC(I)=JN(I,J1)
10 JC(I+3)=JN(I,J2)
RETURN
END


C （七）求单元常数
SUBROUTINE SCL(M,NE,NJ,BL,SI,C0,JE,X,Y)
DIMENSION JE(2,NE),X(NJ),Y(NJ)
REAL BL,SI,C0,DX,DY
J1=JE(1,M)
J2=JE(2,M)
DX=X(J2)-X(J1)
DY=Y(J2)-Y(J1)
BL=SQRT(DX*DX+DY*DY)
SI=DY/BL
C0=DX/BL
RETURN
END

C （八）形成单元刚度矩阵
SUBROUTINE ESM(M,NE,BL,EA,EI,KD)
DIMENSION EA(NE),EI(NE)
REAL KD(6,6),BL,S,G,G1,G2,G3
G=EA(M)/BL
G1=2.0*EI(M)/BL
G2=3.0*G1/BL
G3=2.0*G2/BL
DO 10 I=1,6
DO 10 J=1,6
10 KD(I,J)=0.0
KD(1,1)=G
KD(1,4)=-G
KD(4,4)=G
KD(2,2)=G3
KD(5,5)=G3
KD(2,5)=-G3
KD(2,3)=G2
KD(2,6)=G2
KD(3,5)=-G2
KD(5,6)=-G2
KD(3,3)=2.0*G1
KD(6,6)=2.0*G1
KD(3,6)=G1
DO 20 I=1,5
I1=I+1
DO 30 J=I1,6
30 KD(J,I)=KD(I,J)
20 CONTINUE
RETURN
END

C （九）形成单元坐标转换矩阵
SUBROUTINE CTM(SI,C0,T)
REAL T(6,6),SI,C0
DO 10 I=1,6
DO 10 J=1,6
10 T(I,J)=0.0
T(1,1)=C0
T(1,2)=SI
T(2,1)=-SI
T(2,2)=C0
T(3,3)=1.0
DO 20 I=1,3
DO 20 J=1,3
20 T(I+3,J+3)=T(I,J)
RETURN
END

C （十）形成单元固端力
SUBROUTINE EFX(I,NPF,BL,PF,F0)
DIMENSION PF(4,NPF)
REAL F0(6),A,B,C,G,Q,S,BL
IND=PF(2,I)
A=PF(3,I)
Q=PF(4,I)
C=A/BL
G=C*C
B=BL-A
DO 5 J=1,6
5 F0(J)=0.0
GOTO (10,20,30,40,50,60,70),IND
10 S=Q*A*0.5
F0(2)=-S*(2.0-2.0*G+C*G)
F0(5)=-S*G*(2.0-C)
S=S*A/6.0
F0(3)=-S*(6.0-8.0*C+3.0*G)
F0(6)=S*C*(4.0-3.0*C)
WRITE(*,*) F0(3),F0(6)
GOTO 100
20 S=B/BL
F0(2)=-Q*S*S*(1.0+2.0*C)
F0(5)=-Q*G*(1.0+2.0*S)
F0(3)=-Q*S*S*A
F0(6)=Q*B*G
GOTO 100
30 S=B/BL
F0(2)=6.0*Q*C*S/BL
F0(5)=-F0(2)
F0(3)=Q*S*(2.0-3.0*S)
F0(6)=Q*C*(2.0-3.0*C)
GOTO 100
40 S=Q*A*0.25
F0(2)=-S*(2.0-3.0*G+1.6*G*C)
F0(5)=-S*G*(3.0-1.6*C)
S=S*A
F0(3)=-S*(2.0-3.0*C+1.2*G)/1.5
F0(6)=S*C*(1.0-0.8*C)
GOTO 100
50 F0(1)=-Q*A*(1.0-0.5*C)
F0(4)=-0.5*Q*C*A
GOTO 100
60 F0(1)=-Q*B/BL
F0(4)=-Q*C
GOTO 100
70 S=B/BL
F0(2)=Q*C*(3.0*S+C)
F0(5)=-F0(2)
S=S*B/BL
F0(3)=-Q*S*A
F0(6)=Q*G*B
100 RETURN
END

感谢楼主的提供的资料