[Fortran] 纯文本查看 复制代码
program main
common NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION X(100),Y(100),NCO(200),PROP(100),IB(60),TK(200,200),&
& AL(200),FORC(100),REAC(200),ELST(4,4),V(20)
CHARACTER(len=20) FILE1,FILE2
!初始化程序参数
NRMX=200
NCMX=20
NDF=2
NNE=2
NDFEL=NDF*NNE
IN=5
IO=6
!打开输入输出文件
READ(*,*)FILE1,FILE2
OPEN(UNIT=IO,FILE=FILE1,FORM='FORMATTED',STATUS='UNKNOWN')
OPEN(UNIT=IN,FILE=FILE2,FORM='FORMATTED',STATUS='OLD')
!数据输入
CALL INPUT(X,Y,NCO,PROP,AL,IB,REAC)
!组装总体刚度矩阵
CALL ASSEM(X,Y,NCO,PROP,TK,ELST,AL)
!引入边界条件
CALL BOUND(TK,AL,REAC,IB)
!总体方程求解
CALL SLBSI(TK,AL,V,N,MS,NRMX,NCMX)
!计算单元内力与结点合力
CALL FORCE(NCO,PROP,FORC,REAC,X,Y,AL)
!输出
CALL OUTPT(AL,FORC,REAC)
STOP
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE INPUT(X,Y,NCO,PROP,AL,IB,REAC)
COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION X(1),Y(1),PROP(1),AL(1),IB(1),REAC(1),W(3),IC(2),NCO(1)
WRITE(IO,20)
20 FORMAT('',70('*'))
!读入基本常数
READ(IN,*)NN,NE,NLN,NBN,E
WRITE(IO,21)NN,NE,NLN,NBN,E
21 FORMAT(//'INTERNAL DATA'//'NUMBER OF NODES:',I5/&
& 'NUMBER OF ELEMENTS:',I5/'NUMBER OF LOADED NODES:',&
& I5/'NUMBER OF SUPPORT NODES:',I5/'MODULUS OF ELASTICITY:',&
& F15.0//'NODAL COORDINATES'/7X,'NODE',6X,'X',9X,'Y' )
!输入结点x和y坐标值
READ(IN,*)(I,X(I),Y(I),I=1,NN)
WRITE(IO,2)(I,X(I),Y(I),I=1,NN)
2 FORMAT(I10,2F10.2)
!输入结点号和截面几何性质
WRITE(IO,22)
22 FORMAT(/'ELEMENT CONNECTIVITY AND PROPERTIES'/4X,'ELEMENT',3X,&
& 'START NODE END NODE',5X,'AREA')
DO I=1,NE
READ(IN,*)NUM,IC(1),IC(2),PROP(NUM)
WRITE(IO,34)NUM,IC(1),IC(2),PROP(NUM)
N1=NNE*(NUM-1)
NCO(N1+1)=IC(1)
NCO(N1+2)=IC(2)
END DO
34 FORMAT(3I10,F15.5)
N=NN*NDF
DO I=1,N
REAC(I)=0
AL(I)=0
END DO
WRITE(IO,23)
23 FORMAT(/'NODAL LOADS'/7X,'NODE',5X,'PX',8X,'PY')
DO I=1,NLN
READ(IN,*)NUM,(W(K),K=1,NDF)
WRITE(IO,2)NUM,(W(K),K=1,NDF)
DO K=1,NDF
L=NDF*(NUM-1)+K
AL(L)=W(K)
END DO
END DO
!READ BOUNDARY DATA.STORE UNKNOWN STATUS INDICATORS
WRITE(IO,24)
24 FORMAT(/'BOUNDARY CONDITION DATA'/23X,'STATUS',14X,'PRESCRIBED VALUES'/15X,'(0:PRESCRIBED,1:FREE)'/7X,'NODE',8X,'U',9X,&
& 'V',16X,'U',9X,'V')
DO I=1,NBN
READ(IN,*)NUM,(IC(K),K=1,NDF),(W(K),K=1,NDF)
WRITE(IO,9) NUM,(IC(K),K=1,NDF),(W(K),K=1,NDF)
L1=(NDF+1)*(I-1)+1
L2=NDF*(NUM-1)
IB(L1)=NUM
DO K=1,NDF
N1=L1+K
N2=L2+K
IB(N1)=IC(K)
REAC(N2)=W(K)
END DO
END DO
9 FORMAT(3I10,10X,2F10.4)
RETURN
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE ASSEM(X,Y,NCO,PROP,TK,ELST,AL)
COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION X(1),Y(1),NCO(1),TK(200,20),ELST(NDFEL,NDFEL),PROP(1),AL(1)
!COMPUTE HALF BAND WIDTH AND STORE IN MS
N1=NNE-1
MS=0
DO I=1,NE
L1=NNE*(I-1)
DO J=1,N1
L2=L1+J
J1=J+1
DO K=J1,NNE
L3=L1+K
L=ABS(NCO(L2)-NCO(L3))
IF((MS-L).LT.0)THEN
MS=L
ENDIF
ENDDO
ENDDO
ENDDO
MS=NDF*(MS+1)
!CLEAR THE TOTAL STIFFENSS MATRIX
DO I=1,N
DO J=1,MS
TK(I,J)=0
END DO
END DO
DO NEL=1,NE
CALL STIFF(NEL,X,Y,NCO,PROP,ELST,AL)
CALL ELASS(NEL,NCO,TK,ELST)
END DO
RETURN
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE STIFF(NEL,X,Y,NCO,PROP,ELST,AL)
COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION X(1),Y(1),NCO(1),PROP(1),ELST(4,4),AL(1)
L=NNE*(NEL-1)
N1=NCO(L+1)
N2=NCO(L+2)
D=SQRT((X(N2)-X(N1))**2+(Y(N2)-Y(N1))**2)
CO=(X(N2)-X(N1))/D
SI=(Y(N2)-Y(N1))/D
COFE=E*PROP(NEL)/D
ELST(1,1)=COFE*CO*CO
ELST(1,2)=COFE*CO*SI
ELST(2,2)=COFE*SI*SI
DO I=1,2
DO J=1,2
K1=I+NDF
K2=J+NDF
ELST(K1,K2)=ELST(I,J)
ELST(I,K2)=-ELST(I,J)
ENDDO
ENDDO
ELST(2,3)=-ELST(1,2)
RETURN
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE ELASS(NEL,NCO,TM,ELMAT)
COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION NCO(1),TM(200,20),ELMAT(NDFEL,NDFEL)
L1=NNE*(NEL-1)
DO I=1,NNE
L2=L1+I
N1=NCO(L2)
I1=NDF*(I-1)
J1=NDF*(N1-1)
DO J=1,NNE
L2=L1+J
N2=NCO(L2)
I2=NDF*(J-1)
J2=NDF*(N2-1)
DO K=1,NDF
KI=1
IF((N1-N2).EQ.0)THEN
KI=K
ENDIF
IF((N1-N2).LE.0)THEN
KR=J1+K
IC=J2-KR+1
K1=I1+K
ELSE
KR=J2+K
IC=J1-KR+1
K2=I2+K
ENDIF
DO L=KI,NDF
KC=IC+L
IF((N1-N2).LE.0)THEN
K2=I2+L
ELSE
K1=I1+L
ENDIF
TM(KR,KC)=TM(KR,KC)+ELMAT(K1,K2)
ENDDO
ENDDO
ENDDO
ENDDO
RETURN
END
SUBROUTINE BOUND(TK,AL,REAC,IB)
COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION AL(1),IB(1),REAC(1),TK(200,20)
DO L=1,NBN
L1=(NDF+1)*(L-1)+1
NO=IB(L1)
K1=NDF*(NO-1)
DO I=1,NDF
L2=L1+I
IF(IB(L2).EQ.0)THEN
KR=K1+I
DO J=2,MS
KV=KR+J-1
IF((N-KV).GE.0)THEN
AL(KV)=AL(KV)-TK(KR,J)*REAC(KR)
TK(KR,J)=0
ENDIF
KV=KR-J+1
IF(KV.GT.0)THEN
AL(KV)=AL(KV)-TK(KV,J)*REAC(KR)
TK(KV,J)=0
ENDIF
ENDDO
TK(KR,1)=1
AL(KR)=REAC(KR)
ENDIF
ENDDO
ENDDO
RETURN
END
SUBROUTINE SLBSI(A,B,D,N,MS,NX,MX)
DIMENSION A(NX,MX),B(NX),D(MX)
N1=N-1
DO K=1,N1
C=A(K,1)
K1=K+1
IF((ABS(C)-0.00000001).LE.0)THEN
WRITE(6,2)K
2 FORMAT('*****SINGULARITY IN ROW',I5)
STOP
ELSE
N1=K1+MS-2
L=MIN(N1,N)
DO J=2,MS
D(J)=A(K,J)
ENDDO
DO J=K1,L
A(K,K2)=A(K,K2)/C
ENDDO
B(K)=B(K)/C
DO I=K1,L
K2=I-K1+2
C=D(K2)
DO J=I,L
K2=J-I+1
K3=J-K+1
A(I,K2)=A(I,K2)-C*A(K,K3)
ENDDO
B(I)=B(I)-C*B(K)
ENDDO
ENDIF
ENDDO
IF((ABS(C)-0.0000001).LE.0)THEN
WRITE(6,7)K
7 FORMAT('*****SINGULARITY IN ROW',I5)
STOP
ELSE
B(N)=B(N)/A(N,1)
DO I=1,N1
K=N-1
K1=K+1
N1=K1+MS-2
L=MIN(N1,N)
DO J=K1,L
K2=J-K+1
B(K)=B(K)-A(K,K2)*B(J)
ENDDO
ENDDO
ENDIF
RETURN
END
SUBROUTINE FORCE(NCO,PROP,FORC,REAC,X,Y,AL)
COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION NCO(1),PROP(1),FORC(1),REAC(1),X(1),Y(1),AL(1)
DO I=1,N
REAC(I)=0
ENDDO
DO NEL=1,NE
L=NNE*(NEL-1)
N1=NCO(L+1)
N2=NCO(L+2)
K1=NDF*(N1-1)
K2=NDF*(N2-1)
D=SQRT((X(N2)-X(N1))**2+(Y(N2)-Y(N1))**2)
CO=(X(N2)-X(N1))/D
SI=(Y(N2)-Y(N1))/D
COEF=E*PROP(NEL)/D
FORC(NEL)=COEF*((AL(K2+1)-AL(K1+1))*CO+(AL(K2+2)-AL(K1+2))*SI)
REAC(K1+1)=REAC(K1+1)-FORC(NEL)*CO
REAC(K1+2)=REAC(K1+2)-FORC(NEL)*SI
REAC(K2+1)=REAC(K2+1)+FORC(NEL)*CO
REAC(K2+2)=REAC(K2+2)+FORC(NEL)*SI
ENDDO
RETURN
END
SUBROUTINE OUTPT(AL,FORC,REAC)
COMMON NRMX,NCMX,NDFEL,NN,NE,NLN,NBN,NDF,NNE,N,MS,IN,IO,E,G
DIMENSION AL(1),REAC(1),FORC(1)
WRITE(IO,1)
1 FORMAT(//1X,70('*')//'RESULTS'//'NODAL DISPLACEMENTS'/7X,'NODE' &
& ,11X,'U',14X,'V')
DO I=1,NN
K1=NDF*(I-1)+1
K2=K1+NDF-1
WRITE(IO,2)I,(AL(J),J=K1,K2)
ENDDO
2 FORMAT(I10,6F15.4)
WRITE(IO,3)
3 FORMAT(/'NODAL REACTIONS'/7X,'NODE',10X,'PX',13X,'PY')
DO I=1,NN
K1=NDF*(I-1)+1
K2=K1+NDF-1
WRITE(IO,2)I,(REAC(J),J=K1,K2)
ENDDO
WRITE(IO,4)
4 FORMAT(/'MEMBER FORCES'/6X,'MEMBER AXIAL FORCE')
DO I=1,NE
WRITE(IO,2)I,FORC(I)
ENDDO
WRITE(IO,5)
5 FORMAT(//1X,70('*'))
RETURN
END