有限元编程算例.docx

1、有限元编程算例有限元编程算例(Fortran)本程序通过Fortran语言编写，程序在Intel Parallel Studio XE 2013 with VS2013中成功运行，程序为计算力学（龙述尧等编）一书中的源程序，仅作研究学习使用，省去了敲写的麻烦。源程序为：!Page149COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)OPEN(5,FI

2、LE=DATAIN)!OPEN(6,FILE=DATAOUT,STATUS=NEW)CALL DATAIF(IND.EQ.0)GOTO 10 EO=EO/(1.0-UN*UN) UN=UN/(1.0-UN)10 CALL TOTSTICALL LOADCALL SUPPORCALL SOLVEQ CALL STRESS PAUSE!STOPENDSUBROUTINE DATACOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,

3、3),S(3,6),TKZ(200,20),EKE(6,6),P(200)READ(5,*)NJ,NE,NZ,NDD,NPJ,INDNJ2=NJ*2NPJ1=NPJ+1READ(5,*)EO,UN,GAMA,TEREAD(5,*)(JM(I,J),J=1,3),I=1,NE)READ(5,*)(CJZ(I,J),J = 1,2),I=1,NJ)!Page150READ(5,*)(NZC(I),I=1,NZ)READ(5,*)(PJ(I,J),J=1,2),I=1,NPJ1)WRITE(6,10)(I,(CJZ(I,J),J=1,2),I=1,NJ)10 FORMAT(4X,2HNO,6X,1H

4、X,6X,1HY/(I6,2X,F7.2,F7.2)RETURNENDSUBROUTINE ELEST(MEO,IASK)COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)CM=CJZ(JE,1)-CJZ(IE,1)BM=CJZ(IE,2)-CJZ(JE,2)CJ=CJZ(IE,1

5、)-CJZ(ME,1)BJ=CJZ(ME,2)-CJZ(IE,2)AE=(BJ*CM-BM*CJ)/2.0IF(IASK.LE.1) GOTO 50DO 10 I=1,3 DO 10 J=1,6 B(I,J)=0.010 CONTINUEB(1,1)=-BJ-BMB(1,3)=BJB(1,5)=BMB(2,2)=-CJ-CMB(2,4)=CJB(2,6)=CMB(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)!Page151B(3,6)=B(1,5)DO 20 I=1,3 DO 20 J=1,6 B(I,J)=B(

6、I,J)/(2.0*AE)20 CONTINUED(1,1)=EO/(1.0-UN*UN)D(1,2)=EO*UN/(1.0-UN*UN)D(2,1)=D(1,2)D(2,2)=D(1,1)D(1,3)=0.0D(2,3)=0.0D(3,1)=0.0D(3,2)=0.0D(3,3)=EO/(2.0*(1.0+UN)DO 30 I=1,3 DO 30 J=1,6 S(I,J)=0.0 DO 30 K=1,3 S(I,J)=S(I,J)+D(I,K)*B(K,J)30 CONTINUEIF(IASK.LE.2) GOTO 50DO 40 I=1,6 DO 40 J=1,6 EKE(I,J)=0.0

7、 DO 40 K=1,3!*Exchange B And S* EKE(I,J)=EKE(I,J)+B(K,I)*S(K,J)*AE*TE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE TOTSTICOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)!Page152DO 20 I=1,NJ2DO 20 J=1,NDDT

8、KZ(I,J)=0.020 CONTINUE!*Not Understanded*DO 30 MEO=1,NE CALL ELEST(MEO,3) DO 30 I=1,3 DO 30 II=1,2 LH=2*(I-1)+II LDH=2*(JM(MEO,I)-1)+II DO 30 J=1,3 DO 30 JJ=1,2 L=2*(J-1)+JJ LZ=2*(JM(MEO,J)-1)+JJ LD=LZ-LDH+1 IF(LD.LE.0) GOTO 30 TKZ(LDH,LD)=TKZ(LDH,LD)+EKE(LH,L)30 CONTINUERETURNENDSUBROUTINE LOADCOMM

9、ON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 10 I=1,NJ2 P(I)=0.010 CONTINUEIF(NPJ.EQ.0) GOTO 30DO 20 I=1,NPJ I1=I+1 J=IFIX(PJ(I1,2) P(J)=PJ(I1,1)20 CONTINUE30 IF(GAMA.LE.0.0) GOTO 50!Page153DO 40 ME

10、O=1,NE CALL ELEST(MEO,1) PE=-GAMA*AE*TE/3.0 IE=JM(MEO,1) JE=JM(MEO,2) ME=JM(MEO,3) P(2*IE)=P(2*IE)+PE P(2*JE)=P(2*JE)+PE P(2*ME)=P(2*ME)+PE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE SUPPORCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3)

11、,S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 60 I=1,NZ MZ=NZC(I) TKZ(MZ,1)=1.0 DO 10 J=2,NDD TKZ(MZ,J)=0.0 10 CONTINUE IF(MZ-NDD)20,20,30 20 JO=MZ GOTO 40 30 JO=NDD 40 DO 50 J = 2,JO J1=MZ-J TKZ(J1+1,J)=0.0 50 CONTINUE P(MZ)=0.0 60 CONTINUERETURNEND!Page154SUBROUTINE SOLVEQCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ

12、2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)NJ1=NJ2-1DO 50 K=1,NJ1 IF(NJ2-K-NDD+1)10,10,20 10 IM=NJ2 GOTO 30 20 IM=K+NDD-1 30 K1=K+1 DO 50 I=K1,IM L=I-K+1 C=TKZ(K,L)/TKZ(K,1) LD1=NDD-L+1 DO 40 J=1,LD1 M=J+I-K TKZ(I,J)=TKZ(I,J)-C*

13、TKZ(K,M) 40 CONTINUE P(I)=P(I)-C*P(K)50 CONTINUEP(NJ2)=P(NJ2)/TKZ(NJ2,1)DO 100 I1 = 1,NJ1I=NJ2-I1!*下面一行可能出错IF(NDD-NJ2+I-1)60,60,7060 JO=NDD GOTO 8070 JO=NJ2-I+180 DO 90 J=2,JOLH=J+I-1P(I)=P(I)-TKZ(I,J)*P(LH)90 CONTINUEP(I)=P(I)/TKZ(I,1)100 CONTINUE!Page155WRITE(6,110)(I,P(2*I-1),P(2*I),I=1,NJ)!*110

14、FORMAT(2X,3HJD=,3X,2HU=,12X,2HV=/(I4,3X,F16.7,3X,F16.7)RETURNENDSUBROUTINE STRESSCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)DIMENSION WY(6),YL(3)DO 60 MEO=1,NECALL ELEST(MEO,2)DO 10 I=1,3DO 10 J=1,2LH=2*(I-1)+JLDH=2*(JM(MEO,I)-1)+JWY(LH)=P(LDH)10