python有限元后处理 fortran有限元

转载

mob6454cc6172e5 2023-07-27 00:52:47

文章标签 python有限元后处理 fortran有限元 P4 子程序程序段 文章分类 Python 后端开发

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*SEK2(3,1)=-S*P
 EK2(3,2)=C*P
 EK2(3,3)=2.0*QEK2(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

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。