有限元编程算例(fortran).docx

_有限元编程算例(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,FILE=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),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,1HX,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)-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(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 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,NDDTKZ(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 LOADCOMMON/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 MEO=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),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,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)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*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 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 CONTINUEDO 20 I=1,3YL(I)=0.0DO 20 J=1,6YL(I)=YL(I)+S(I,J)*WY(J)20 CONTINUESIGX=YL(1)SIGY=YL(2)TOXY=YL(3)PYL=(SIGX+SIGY)/2.0SIG=(SIGX-SIGY)*2/4.0+TOXY*TOXYRYL=SQRT(SIG)SIG1=PYL+RYLSIG2=PYL-RYLIF(SIGY.EQ.SIG2) GOTO 30CETA1=TOXY/(SIGY-SIG2)CETA=90.0-57.29578*ATAN(CETA1)GOTO 40!Page15630 CETA=0.040 WRITE(6,50)MEO,SIGX,SIGY,TOXY,SIG1,SIG2,CETA50 FORMAT(4X,2HE=,I3/2X,3HSX=,F11.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HS1=,F11.3,3X,3HS2=,F11.3,3X,4HCET=,F11.3)!50 FORMAT(4X,2HE=,I3/2X,3HSX=,Fll.3,3X,3HSY=,F11.3,3X,4HTAU=,F11.3/2X,3HSl=,Fll.3,3X,3HS2=,F11.3,3X,4HCET=,F11.3)60 CONTINUE RETURNEND输入文件为datain28,36,9,10,4,01,0.17,0,11,5,22,5,62,6,33,6,73,7,44,7,85,9,66,9,106,10,77,10,117,11,88,11,129,13,1010,13,1410,14,1111,14,1511,15,1212,15,1613,17,1414,17,1814,18,1515,18,1915,19,1616,19,2017,21,1818,21,2218,22,1919,22,2319,23,2020,23,2421,25,2222,25,2622,26,2323,26,2723,27,2424,27,280,61,62,63,60,51,52,53,50,41,42,43,40,31,32,33,30,21,22,23,20,11,12,13,10,01,02,03,07,15,23,31,39,47,49,50,550,0-5E4,2-10E4,4-10E4,6-5E4,8输出结果为：DATAOUTNO X Y 1 0.00 6.00 2 1.00 6.00 3 2.00 6.00 4 3.00 6.00 5 0.00 5.00 6 1.00 5.00 7 2.00 5.00 8 3.00 5.00 9 0.00 4.00 10 1.00 4.00 11 2.00 4.00 12 3.00 4.00 13 0.00 3.00 14 1.00 3.00 15 2.00 3.00 16 3.00 3.00 17 0.00 2.00 18 1.00 2.00 19 2.00 2.00 20 3.00 2.00 21 0.00 1.00 22 1.00 1.00 23 2.00 1.00 24 3.00 1.00 25 0.00 0.00 26 1.00 0.00 27 2.00 0.00 28 3.00 0.00 JD= U= V= 1 -29766.873 -1173917.750 2 -14003.185 -1174018.875 3 -3753.270 -1179518.125 4 0.000 -1181719.750 5 -26382.471 -1072681.500 6 -10746.993 -1073615.000 7 -2064.593 -1082360.750 8 0.000 -1085873.250 9 -13536.995 -964010.125 10 3372.794 -970055.125 11 7268.415 -989269.125 12 0.000 -998401.812 13 7816.581 -835383.438 14 27176.234 -861713.938 15 22063.230 -905726.125 16 0.000 -927165.188 17 29514.479 -665602.875 18 53419.637 -747340.438 19 34876.832 -839806.812 20 0.000 -881219.125 21 29580.273 -416288.719 22 52944.918 -632601.125 23 17504.195 -803765.688 24 0.000 -859481.938 25 0.000 0.000 26 -120102.820 -583505.375 27 -76202.375 -787347.188 28 0.000 -829170.812 E= 1 SX= -1489.530 SY=-101489.383 TAU= -1489.531 S1= -1467.348 S2=-101511.562 CET= 179.147 E= 2 SX= -1475.844 SY=-100654.875 TAU= -1790.500 S1= -1443.531 S2=-100687.188 CET= 178.966 E= 3 SX= -7021.670 SY=-101597.672 TAU= -3741.688 S1= -6873.875 S2=-101745.469 CET= 177.738 E= 4 SX= -8067.500 SY= -98528.750 TAU= -4459.156 S1= -7848.227 S2= -98748.023 CET= 177.185 E= 5 SX= -13143.328 SY= -99391.750 TAU= -1662.500 S1= -13111.293 S2= -99423.781 CET= 178.896 E= 6 SX= -14652.781 SY= -98337.500 TAU= -1501.062 S1= -14625.867 S2= -98364.414 CET= 178.973 E= 7 SX= -2923.122 SY=-109168.297 TAU= -5888.469 S1= -2597.762 S2=-109493.656 CET= 176.837 E= 8 SX= -716.078 SY=-103681.562 TAU= -8617.406 S1= 0.148 S2=-104397.789 CET= 175.249 E= 9 SX= -9188.316 SY=-105121.867 TAU= -9771.594 S1= -8203.125 S2=-106107.062 CET= 174.243 E= 10 SX= -12285.000 SY= -95180.250 TAU= -12199.594 S1= -10526.887 S2= -96938.359 CET= 171.799 E= 11 SX= -14170.516 SY= -95500.750 TAU= -5489.531 S1= -13801.664 S2= -95869.602 CET= 176.156 E= 12 SX= -22797.406 SY= -91347.000 TAU= -3902.844 S1= -22575.914 S2= -91568.492 CET= 176.752 E= 13 SX= -5104.269 SY=-129494.438 TAU= -11708.750 S1= -4011.727 S2=-130586.977 CET= 174.669 E= 14 SX= 969.672 SY=-108176.375 TAU= -21424.750 S1= 5024.582 S2=-112231.281 CET= 169.283 E= 15 SX= -14954.572 SY=-110883.469 TAU= -18383.531 S1= -11552.273 S2=-114285.766 CET= 169.515 E= 16 SX= -19890.141 SY= -86924.312 TAU= -25131.188 S1= -11514.844 S2= -95299.609 CET= 161.569 E= 17 SX= -22109.688 SY= -87301.625 TAU= -10225.406 S1= -20543.453 S2= -88867.859 CET= 171.292 E= 18 SX= -35190.453 SY= -77219.000 TAU= -9162.000 S1= -33280.023 S2= -79129.430 CET= 168.222 E= 19 SX= -9785.850 SY=-171444.172 TAU= -20524.969 S1= -7220.594 S2=-174009.422 CET= 172.876 E= 20 SX= 4594.438 SY=-113592.375 TAU= -46145.688 S1= 20477.398 S2=-129475.336 CET= 161.007 E= 21 SX= -25287.307 SY=-118672.312 TAU= -30023.750 S1= -16467.512 S2=-127492.109 CET= 163.629 E= 22 SX= -30634.422 SY= -71127.188 TAU= -44991.469 S1= -1543.715 S2=-100217.891 CET= 147.114 E= 23 SX= -34259.609 SY= -71743.438 TAU= -14637.906 S1= -29220.699 S2= -76782.344 CET= 161.005 E= 24 SX= -43958.047 SY= -53418.938 TAU= -17697.562 S1= -30369.627 S2= -67007.359 CET= 142.482 E= 25 SX= -19028.160 SY=-252549.000 TAU= -34958.688 S1= -13907.055 S2=-257670.094 CET= 171.666 E= 26 SX= 3973.812 SY=-114063.750 TAU= -92238.344 S1= 54459.047 S2=-164548.984 CET= 151.307 E= 27 SX= -39180.809 SY=-121400.055 TAU= -39312.688 S1= -23409.074 S2=-137171.781 CET= 158.140 E= 28 SX= -42804.766 SY= -43317.938 TAU= -65723.062 S1= 22662.211 S2=-108784.914 CET= 135.112 E= 29 SX= -42224.094 SY= -43219.188 TAU= -10273.375 S1= -32436.225 S2= -53007.055 CET= 136.386 E= 30 SX= -21830.422 SY= -25448.312 TAU= -23810.344 S1= 239.594 S2= -47518.328 CET= 137.172 E= 31 SX= -48815.199 SY=-424587.344 TAU= -79800.078 S1= -32570.844 S2=-440831.688 CET= 168.494 E= 32 SX=-132271.750 SY= -71582.000 TAU=-175409.250 S1= 76087.781 S2=-279941.531 CET= 130.093 E= 33 SX= -45090.102 SY= -56761.105 TAU= 804.781 S1= -45034.867 S2= -56816.336 CET= 3.926 E= 34 SX= 42332.711 SY= -9221.938 TAU= -47066.344 S1= 70218.328 S2= -37107.555 CET= 149.354 E= 35 SX= -20899.344 SY= -19971.375 TAU= 16235.219 S1= -4193.512 S2= -36677.207 CET= 45.819 E= 36 SX= 73163.914 SY= -17873.250 TAU= -17873.344 S1= 76547.250 S2= -21256.586 CET= 169.281Welcome ToDownload !欢迎您的下载，资料仅供参考！精品资料