FORTRAN 开发平面薄壁结构有限元程序



PROGRAM MAIN
CHARACTER*30 INFILE
REAL K
DIMENSION NOPNT(20),NOFIX(20),INFOC(20,3),JAD(16),X(100), PLOAD(20,3),PRESC(20,3),EK1(4,4),EK2(6,6),EK4(8,8)
COMMON/LIMT/ME1,ME2,ME4,MJ
COMMON/CTRL/NE1,NE2,NE4,NJE1,NJE2,NJE4,NJ,NCJ,NPJ,NFJ
COMMON/TOPL/IA1(20,2),IA2(20,2),IA4(20,4),XY(40,2)
COMMON/STIF/K(100,101)
COMMON/C/PROPS(5,3)
COMMON/C1/MATNO1(20)/C2/MATNO2(20)/C4/MATNO4(20),AREA,R(8)
ME1=20
ME2=20
ME4=20
MJ=40
NJE1=2
NJE2=2
NJE4=4
NFJ=3
NPROP=3
WRITE(*,*)'PLEASE ENTER DATA FILE NAME:'
!READ(*,'(A)')INFILE

!OPEN(1,FILE=INFILE,STATUS='OLD')
OPEN(1,FILE="1.txt",STATUS='OLD')
READ(1,*)NE1,NE2,NE4,NJ
READ(1,*)((IA1(I,J),J=1,NJE1),I=1,NE1)
READ(1,*)((IA2(I,J),J=1,NJE2),I=1,NE2)
READ(1,*)((IA4(I,J),J=1,NJE4),I=1,NE4)
READ(1,*)(XY(I,1),XY(I,2),I=1,NJ)
READ(1,*)NMATS
READ(1,*)((PROPS(I,J),J=1,NPROP),I=1,NMATS)
READ(1,*)(MATNO1(I),I=1,NE1)
READ(1,*)(MATNO2(I),I=1,NE2)
READ(1,*)(MATNO4(I),I=1,NE4)
READ(1,*)NCJ
READ(1,*)(NOFIX(I),I=1,NCJ)
READ(1,*)((INFOC(I,J),J=1,NFJ),I=1,NCJ)
READ(1,*)((PRESC(I,J),J=1,NFJ),I=1,NCJ)
READ(1,*)NPJ
READ(1,*)(NOPNT(I),I=1,NPJ)
READ(1,*)((PLOAD(I,J),J=1,NFJ),I=1,NPJ)

write(*,*)NE1,NE2,NE4,NJ
write(*,300)((IA1(I,J),J=1,NJE1),I=1,NE1)
write(*,300)((IA2(I,J),J=1,NJE2),I=1,NE2)
write(*,400)((IA4(I,J),J=1,NJE4),I=1,NE4)
write(*,500)(XY(I,1),XY(I,2),I=1,NJ)
write(*,*)NMATS
write(*,600)((PROPS(I,J),J=1,NPROP),I=1,NMATS)
write(*,*)(MATNO1(I),I=1,NE1)
write(*,*)(MATNO2(I),I=1,NE2)
write(*,*)(MATNO4(I),I=1,NE4)
write(*,*)NCJ
write(*,*)(NOFIX(I),I=1,NCJ)
write(*,*)((INFOC(I,J),J=1,NFJ),I=1,NCJ)
write(*,*)((PRESC(I,J),J=1,NFJ),I=1,NCJ)
write(*,*)NPJ
write(*,*)(NOPNT(I),I=1,NPJ)
write(*,*)((PLOAD(I,J),J=1,NFJ),I=1,NPJ)
write(*,*)"******************************"
300 format(1x,2i5)

400 format(1x,4i5)
500 format(1x,2f20.2)

600 format(1x,3f10.3)
700 format(1x,2f20.2)
800 format(1x,2f20.2)

!goto 1100

CLOSE(1)
NTF=NJ*NFJ
DO II=1,NTF
X(II)=0.0
DO JJ=1,NTF+1
K(II,JJ)=0.0
END DO
END DO
DO IE=1,NE1
CALL ESTIF1(IE,EK1)
CALL EAD(IE,ME1,NFJ,NJE1,2,JAD,IA1)
CALL FORK(4,JAD,EK1)
END DO
DO IE=1,NE2
CALL ESTIF2(IE,EK2)
CALL EAD(IE,ME2,NFJ,NJE2,3,JAD,IA2)
CALL FORK(6,JAD,EK2)
END DO
DO IE=1,NE4
CALL ESTIF4(IE,EK4)
CALL EAD(IE,ME4,NFJ,NJE4,2,JAD,IA4)
CALL FORK(8,JAD,EK4)
END DO
DO I=1,NPJ
DO J=1,3
II=(NOPNT(I)-1)*3+J
K(II,NTF+1)=PLOAD(I,J)
END DO
END DO
DO I=1,NCJ
DO J=1,3
II=(NOFIX(I)-1)*3+J
JJ=INFOC(I,J)
IF(JJ.NE.0)THEN
K(II,II)=K(II,II)*10E18
K(II,NTF+1)=K(II,II)*PRESC(I,J)
END IF
END DO
END DO
K(15,15)=10E18
K(24,24)=10E18
CALL GS(NTF,X)
CALL RESULT(X)
write(*,900)((EK2(I,J),J=1,6),I=1,6)
900 format(1x,6f15.5)
1100   END

!数据输出程序段
SUBROUTINE RESULT(X)
CHARACTER*30 OUTFILE
DIMENSION X(100),EK1(4,4),EK2(6,6),EK4(8,8)
DIMENSION D1(4),D2(6),D4(8),P1(4),P2(6),P4(8)
COMMON/CTRL/NE1,NE2,NE4,NJE1,NJE2,NJE4,NJ,NCJ,NPJ,NFJ
COMMON/TOPL/IA1(20,2),IA2(20,2),IA4(20,4),XY(40,2)
COMMON/C/PROPS(5,3)/C4/MATNO4(20),AREA,R(8)
WRITE(*,*)'ENTER THE OUTFILE NAME:'
!READ(*,'(A)')OUTFILE
!OPEN(2,FILE=OUTFILE,STATUS='UNKNOWN')
OPEN(2,FILE="2.txt",STATUS='UNKNOWN')
WRITE(2,*)'NODAL DISPLACEMENTS'
WRITE(2,*)' UX UY UXY'
DO I=1,NJ
WRITE(2,'(1X,I3,2X,3F11.5)')I,(X((I-1)*3+J),J=1,NFJ)
END DO
WRITE(2,*)'NODAL FORCES ON BAR ELEMENTS:'
WRITE(2,*)'PXI PYI PXJ PYJ'
DO IE=1,NE1
CALL ESTIF1(IE,EK1)
I=IA1(IE,1)
D1(1)=X((I-1)*3+1)
D1(2)=X((I-1)*3+2)
I=IA1(IE,2)
D1(3)=X((I-1)*3+1)
D1(4)=X((I-1)*3+2)
DO I=1,4
P1(I)=0.
 DO J=1,4
P1(I)=P1(I)+EK1(I,J)*D1(J)
END DO
END DO
WRITE(2,'(1X,I3,4F11.3)')IE,(P1(I),I=1,4)
END DO
WRITE(2,*)'NODAL FORCES ON BIAM ELEMENTS:'
WRITE(2,1)
1 FORMAT(12X,'PXI',8X,'PYI',8X,'MI',9X,'PXJ',8X,'PYJ',8X,'MJ')
DO IE=1,NE2
CALL ESTIF2(IE,EK2)
I=IA2(IE,1)
D2(1)=X((I-1)*3+1)
D2(2)=X((I-1)*3+2)
D2(3)=X((I-1)*3+3)
I=IA2(IE,2)
D2(4)=X((I-1)*3+1)
D2(5)=X((I-1)*3+2)
D2(6)=X((I-1)*3+3)
DO I=1,6
P2(I)=0.
DO J=1,6
P2(I)=P2(I)+EK2(I,J)*D2(J)
END DO
END DO
WRITE(2,'(1X,I3,6F11.3)')IE,(P2(I),I=1,6)
END DO
WRITE(2,*)'SHEAR STRESS AND NODAL FORCES ON SHEAR-ONLY PLATE ELEMENTS:'
WRITE(2,*)
2 FORMAT(9X,'STR',7X,'PXI',7X,'PYI',7X,'PXJ',7X,'PYJ',7X,'PXM' ,7X,'PYM',7X,'PXN',7X,'PYN')
DO IE=1,NE4
CALL ESTIF4(IE,EK4)
I=IA4(IE,1)
D4(1)=X((I-1)*3+1)
D4(2)=X((I-1)*3+2)
I=IA4(IE,2)
D4(3)=X((I-1)*3+1)
D4(4)=X((I-1)*3+2)
I=IA4(IE,3)
D4(5)=X((I-1)*3+1)
D4(6)=X((I-1)*3+2)
I=IA4(IE,4)
D4(7)=X((I-1)*3+1)
D4(8)=X((I-1)*3+2)
LMAT=MATNO4(IE)
G4=PROPS(LMAT,1)
TAU=0.
DO I=1,8
TAU=TAU+R(I)*D4(I)
END DO
TAU=TAU*0.5*G4/AREA
DO I=1,8
P4(I)=0.
DO J=1,8
P4(I)=P4(I)+EK4(I,J)*D4(J)
END DO
END DO
WRITE(2,'(1X,I3,F9.3,8F10.3)')IE,TAU,(P4(I),I=1,8)
END DO
CLOSE(2)

END

!梁无刚度矩阵子程序
 
subroutine ESTIF2(MM,EK2)
dimension EK2(6,6)
common/TOPl/IA1(20,2),IA2(20,2),IA4(20,4),XY(40,2)
common/C/PROPS(5,3)/C2/MATNO2(20)
LMAT=MATNO2(MM)
E2=PROPS(LMAT,1)
F2=PROPS(LMAT,2)
AJ2=PROPS(LMAT,3)
I=IA2(MM,1)
X1=XY(I,1)
Y1=XY(I,2)
I=IA2(MM,2)
X2=XY(I,1)
Y2=XY(I,2)
B=sqrt((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1))
C=(X2-X1)/B
S=(Y2-Y1)/B
F=12.0*E2*AJ2/(B*B*B)
H=E2*F2/B
P=6.0*E2*AJ2/(B*B)
Q=2.0*E2*AJ2/B

EK2(1,1)=F*S*S+H*C*C
EK2(2,1)=S*C*(H-F)
EK2(2,2)=F*C*C+H*S*S

EK2(3,1)=-S*P
EK2(3,2)=C*P
EK2(3,3)=2.0*Q

EK2(4,1)=-EK2(1,1)
EK2(4,2)=-EK2(2,1)
EK2(4,3)=S*P
EK2(4,4)=EK2(1,1)

EK2(4,1)=-EK2(1,1)
EK2(4,2)=-EK2(2,1)
EK2(4,3)=S*P
EK2(4,4)=EK2(1,1)

EK2(5,1)=-EK2(2,1)
EK2(5,2)=-EK2(2,2)
EK2(5,3)=-C*P
EK2(5,4)=EK2(2,1)
EK2(5,5)=EK2(2,2)

EK2(6,1)=-S*P
EK2(6,2)=C*P
EK2(6,3)=Q
EK2(6,4)=S*P
EK2(6,5)=-C*P
EK2(6,6)=2.0*Q
DO 1 I=2,6
 DO 1 J=1,I-1
1 EK2(J,I)=EK2(I,J)
END

!单元地址程序段
SUBROUTINE EAD(IE,ME,NFJ,NJE,NFJE,JAD,IA)
DIMENSION JAD(16),IA(ME,NJE)
DO I=1,NJE
DO J=1,NFJE
IJ=NFJE*(I-1)+J
II=NFJ*(IA(IE,I)-1)+J
JAD(IJ)=II
END DO
END DO
END

!总体刚度矩阵程序段
SUBROUTINE FORK(NM,JAD,EK)
REAL K
DIMENSION EK(NM,NM),JAD(16)
COMMON/STIF/K(100,101)
DO I=1,NM
DO J=1,NM
II=JAD(I)
JJ=JAD(J)
K(II,JJ)=K(II,JJ)+EK(I,J)
END DO
END DO
END
!高斯消去法子程序
SUBROUTINE GS(N,X)
DIMENSION X(N)
COMMON/STIF/A(100,101)
DO 10 K=1,N-1
DO 10 I=K+1,N
DO 10 J=K+1,N+1
10 A(I,J)=A(I,J)-A(I,K)*A(K,J)/A(K,K)
X(N)=A(N,N+1)/A(N,N)
DO 30 K=N-1,1,-1
X(K)=0.0
DO 20 J=K+1,N
20 X(K)=X(K)+A(K,J)*X(J)
30 X(K)=(A(K,N+1)-X(K))/A(K,K)
END

!杆单元刚度矩阵子程序
SUBROUTINE ESTIF1(MM,EK1)
REAL LIJ
DIMENSION EK1(4,4)
COMMON/TOPL/IA1(20,2),IA2(20,2),IA4(20,4),XY(40,2)
COMMON/C/PROPS(5,3)/C1/MATNO1(20)
LMAT=MATNO1(MM)
E1=PROPS(LMAT,1)
F1=PROPS(LMAT,2)
I=IA1(MM,1)
X1=XY(I,1)
Y1=XY(I,2)
I=IA1(MM,2)
X2=XY(I,1)
Y2=XY(I,2)
LIJ=SQRT((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1))
C=(X2-X1)/LIJ
S=(Y2-Y1)/LIJ
F=F1*E1/LIJ
EK1(1,1)=F*C*C
EK1(2,1)=F*C*S
EK1(2,2)=F*S*S
EK1(3,1)=-F*C*C
EK1(3,2)=-F*C*S
EK1(3,3)=F*C*C
EK1(4,1)=-F*C*S
EK1(4,2)=-F*S*S
EK1(4,3)=F*S*C
EK1(4,4)=F*S*S
DO 1 I=1,3
DO 1 J=I+1,4
1 EK1(I,J)=EK1(J,I)
END

SUBROUTINE ESTIF4(MM,EK4)
REAL L12
DIMENSION EK4(8,8),X(4),Y(4)
COMMON/TOPL/IA1(20,2),IA2(20,2),IA4(20,4),XY(40,2)
COMMON/C/PROPS(5,3)/C4/MATNO4(20),AREA,R(8)
LMAT=MATNO4(MM)
G4=PROPS(LMAT,1)
T4=PROPS(LMAT,2)
DO I=1,4
II=IA4(MM,I)
X(I)=XY(II,1)
Y(I)=XY(II,2)
END DO
A1=X(2)*Y(3)+Y(2)*X(1)-X(2)*Y(1)-X(3)*Y(2)
A2=X(3)*Y(4)+X(4)*Y(1)-X(4)*Y(3)-X(1)*Y(4)
AREA=0.5*(A1+A2)
AA=G4*T4/(4.*AREA)
X21=X(2)-X(1)
Y21=Y(2)-Y(1)
X41=X(4)-X(1)
Y41=Y(4)-Y(1)
X42=X(4)-X(2)
Y42=Y(4)-Y(2)
X31=X(3)-X(1)
Y31=Y(3)-Y(1)
L12=SQRT(X21*X21+Y21*Y21)
C=X21/L12
S=Y21/L12
R(1)=-S*S*X41+C*C*X42+C*S*Y41+C*S*Y42
R(2)=C*S*X41+C*S*X42-C*C*Y41+S*S*Y42
R(3)=(S*S-C*C)*X31-2*S*C*Y31
R(4)=-2*C*S*X31+(C*C-S*S)*Y31
DO I=5,8
R(I)=-R(I-4)
END DO
DO I=1,8
DO J=1,8
EK4(I,J)=R(I)*R(J)*AA
END DO
END DO
END