1424.f90

有限元分析程序

MODULE ALLDATA
  PARAMETER (MAXNODE=500)
  PARAMETER (MAXELEM=500)
  PARAMETER (MAXFORC=500)
  PARAMETER (MAXMATO=3000)
  INTEGER:: NJ,NM,NE,NF
  INTEGER:: NR(MAXNODE,6)
  REAL*8 :: X(MAXNODE),Y(MAXNODE),Z(MAXNODE)
  REAL*8 :: XE(MAXELEM,7)
  INTEGER:: JE(MAXELEM,5)
  INTEGER:: NJF(MAXFORC),NDF(MAXFORC)
  REAL*8 :: VF(MAXFORC)
  INTEGER:: ID(MAXNODE,6),NEQ
  REAL*8 :: AK(MAXMATO,MAXMATO)
END MODULE ALLDATA

PROGRAM MAIN
USE ALLDATA      
CHARACTER*12 NAME0,NAME1
WRITE(*,*) 'INPUT FILE NAME ---'      
READ*,NAME0
WRITE(*,*) 'OUTPUT FILE NAME ---'      
READ*,NAME1
OPEN(5,FILE=NAME0)
OPEN(6,FILE=NAME1)
CALL INPUT
CALL FORMID
CALL FORMK
CALL DISPL
STOP
END

SUBROUTINE INPUT
USE ALLDATA
READ(5,*) NJ,NM,NE,NF

DO I=1,NJ
READ(5,*) II,(NR(I,J),J=1,6),X(I),Y(I),Z(I)
ENDDO
DO I=1,NM
READ(5,*) II,(XE(I,J),J=1,7)
ENDDO
DO I=1,NE
READ(5,*) II,(JE(I,J),J=1,4)
ENDDO
DO I=1,NF
READ(5,*) II,NJF(I),NDF(I),VF(I)
ENDDO

WRITE(6,1000)
WRITE(6,1010)
WRITE(6,1020)
WRITE(6,1030)
WRITE(6,1040)
WRITE(6,*)
WRITE(6,1050) NJ
DO I=1,NJ
WRITE(6,1060) I,(NR(I,J),J=1,6),X(I),Y(I),Z(I)
ENDDO
WRITE(6,*)
WRITE(6,1070) NM
DO I=1,NM
WRITE(6,1080) I,(XE(I,J),J=1,7)
ENDDO
WRITE(6,*)
WRITE(6,1090) NE
DO I=1,NE
WRITE(6,1100) I,(JE(I,J),J=1,4)
ENDDO
WRITE(6,*)
WRITE(6,1110) NF
DO I=1,NF
WRITE(6,1120) I,NJF(I),NDF(I),VF(I)
ENDDO

1000  FORMAT(20X,'EEEEE     SSSS      A      PPPP ')
1010  FORMAT(20X,'E        S         A A     P   P')
1020  FORMAT(20X,'EEEEE    SSSSS    AAAAA    PPPP ')
1030  FORMAT(20X,'E            S    A   A    P    ')
1040  FORMAT(20X,'EEEEE    SSSS     A   A    P    ')
1050  FORMAT('NUMBER OF NODE=',I5)
1060  FORMAT(7I5,3E13.5)
1070  FORMAT('NUMBER OF MATERAL=',I5)
1080  FORMAT(I5,7E13.5)
1090  FORMAT('NUMBER OF ELEMENT=',I5)
1100  FORMAT(6I5)
1110  FORMAT('NUMBER OF FORCE=',I5)
1120  FORMAT(3I5,E13.5)

RETURN
END

SUBROUTINE FORMID
USE ALLDATA

K=0
DO I=1,NJ
DO J=1,6
IF (NR(I,J)==0) THEN
K=K+1
ID(I,J)=K
ENDIF
ENDDO
ENDDO
NEQ=K

WRITE(6,*)
WRITE(6,1000)
DO I=1,NJ
WRITE(6,1010) I,(ID(I,J),J=1,6)
ENDDO

1000  FORMAT('ORDER FOR DEGREE OF FREEDOM')
1010  FORMAT(7I6)

RETURN
END

SUBROUTINE FORMK
USE ALLDATA
REAL*8 ::SK(12,12),SM(12,12)
INTEGER::M(12)

AK(1:NEQ,1:NEQ)=0

DO I=1,NE
CALL LY(I,SK,SM)
M(1) =ID(JE(I,1),1)
M(2) =ID(JE(I,1),2)
M(3) =ID(JE(I,1),3)
M(4) =ID(JE(I,1),4)
M(5) =ID(JE(I,1),5)
M(6) =ID(JE(I,1),6)
M(7) =ID(JE(I,2),1)
M(8) =ID(JE(I,2),2)
M(9) =ID(JE(I,2),3)
M(10)=ID(JE(I,2),4)
M(11)=ID(JE(I,2),5)
M(12)=ID(JE(I,2),6)
DO J=1,12
DO K=1,12
IF (M(J)*M(K)>0) AK(M(J),M(K))=AK(M(J),M(K))+SK(J,K)
ENDDO
ENDDO
ENDDO

WRITE(6,*)
WRITE(6,1000)
WRITE(6,1010) NEQ
DO I=1,NEQ
DO J=1,I
IF (ABS(AK(I,J))>0.000001) THEN
WRITE(6,1020) I,J,AK(I,J)
ENDIF
ENDDO
ENDDO

1000  FORMAT('WHOLLY MASS MATRIX AND STIFFNESS MATRIX')
1010  FORMAT('NUMBER OF EQUATIONS=',I5)
1020  FORMAT(2I5,E18.8)

RETURN
END

SUBROUTINE LY(IM,SK,SM)
USE ALLDATA
USE numerical_libraries
    
REAL*8 ::SK(12,12),SM(12,12),SKE(12,12),SME(12,12),R(12,12)
REAL*8 ::L12X,L12Y,L12Z,L13X,L13Y,L13Z
      
II=JE(IM,1)
JJ=JE(IM,2)
KK=JE(IM,3)
DX=X(JJ)-X(II)
DY=Y(JJ)-Y(II)
DZ=Z(JJ)-Z(II)
DL=SQRT(DX**2+DY**2+DZ**2)

E  =XE(JE(IM,4),1)
G  =XE(JE(IM,4),2)
RO =XE(JE(IM,4),3)
A  =XE(JE(IM,4),4)
DIX=XE(JE(IM,4),5)
DIY=XE(JE(IM,4),6)
DIZ=XE(JE(IM,4),7)

DO I=1,12
DO J=1,12
SKE(I,J)=0.0
SME(I,J)=0.0
SK(I,J)=0.0
SM(I,J)=0.0
R(I,J)=0.0
ENDDO
ENDDO

SKE( 1, 1)=E*A/DL
SKE( 2, 2)=12*E*DIZ/DL/DL/DL
SKE( 3, 3)=12*E*DIY/DL/DL/DL
SKE( 4, 4)=G*DIX/DL
SKE( 5, 3)=-6*E*DIY/DL/DL
SKE( 5, 5)=4*E*DIY/DL
SKE( 6, 2)=6*E*DIZ/DL/DL
SKE( 6, 6)=4*E*DIZ/DL
SKE( 7, 1)=-E*A/DL
SKE( 7, 7)=E*A/DL
SKE( 8, 2)=-12*E*DIZ/DL/DL/DL
SKE( 8, 6)=-6*E*DIZ/DL/DL
SKE( 8, 8)=12*E*DIZ/DL/DL/DL
SKE( 9, 3)=-12*E*DIY/DL/DL/DL
SKE( 9, 5)=6*E*DIY/DL/DL
SKE( 9, 9)=12*E*DIY/DL/DL/DL
SKE(10, 4)=-G*DIX/DL
SKE(10,10)=G*DIX/DL
SKE(11, 3)=-6*E*DIY/DL/DL
SKE(11, 5)=2*E*DIY/DL
SKE(11, 9)=6*E*DIY/DL/DL
SKE(11,11)=4*E*DIY/DL
SKE(12, 2)=6*E*DIZ/DL/DL
SKE(12, 6)=2*E*DIZ/DL
SKE(12, 8)=-6*E*DIZ/DL/DL
SKE(12,12)=4*E*DIZ/DL
DO I=1,11
DO J=I+1,12
SKE(I,J)=SKE(J,I)
ENDDO
ENDDO

R(1,1)=DX/DL
R(1,2)=DY/DL
R(1,3)=DZ/DL
G1=X(KK)-X(II)
G2=Y(KK)-Y(II)
G3=Z(KK)-Z(II)
GL=SQRT(G1*G1+G2*G2+G3*G3)
G1=G1/GL
G2=G2/GL
G3=G3/GL
A1=R(1,2)*G3-R(1,3)*G2
A2=R(1,3)*G1-R(1,1)*G3
A3=R(1,1)*G2-R(1,2)*G1
A=SQRT(A1*A1+A2*A2+A3*A3)
R(3,1)=A1/A
R(3,2)=A2/A
R(3,3)=A3/A
R(2,1)=(R(1,3)*A2-R(1,2)*A3)/A
R(2,2)=(R(1,1)*A3-R(1,3)*A1)/A
R(2,3)=(R(1,2)*A1-R(1,1)*A2)/A

DO I=1,3
R(3*I+1,3*I+1)=R(1,1)
R(3*I+1,3*I+2)=R(1,2)
R(3*I+1,3*I+3)=R(1,3)
R(3*I+2,3*I+1)=R(2,1)
R(3*I+2,3*I+2)=R(2,2)
R(3*I+2,3*I+3)=R(2,3)
R(3*I+3,3*I+1)=R(3,1)
R(3*I+3,3*I+2)=R(3,2)
R(3*I+3,3*I+3)=R(3,3)
ENDDO

SK=MATMUL(TRANSPOSE(R),MATMUL(SKE,R))
RETURN
END

SUBROUTINE DISPL
USE ALLDATA
USE numerical_libraries

REAL*8 ::F(MAXMATO),DIS(MAXMATO),DISN(MAXNODE,6)

DO I=1,NF
F(ID(NJF(I),NDF(I)))=VF(I)
ENDDO

CALL DLSARG(NEQ,AK,MAXMATO,F,1,DIS)

DO I=1,NJ
DO J=1,6
DISN(I,J)=0.0
IF(ID(I,J)>0) DISN(I,J)=DIS(ID(I,J))
ENDDO
ENDDO

WRITE(6,*)
WRITE(6,1000)
WRITE(6,1010)
DO I=1,NJ
WRITE(6,1020) I,(DISN(I,J),J=1,6)
ENDDO

1000  FORMAT('NODE DISPLACEMENT RESULT:')
1010  FORMAT(' NODE ',9X,'X',14X,'Y',13X,'Z',14X,'RX',13X,'RY',13X,'RZ')
1020  FORMAT(I4,2X,6E15.6)

RETURN
END