有限元编程算例

1、有限元编程算例(Fortran)本程序通过 Fortran 语言编写，程序在 Intel Parallel Studio XE 2013 withVS2013中成功运行，程序为计算力学(龙述尧等编)一书中的源程序，仅作研究学习使用，省去了敲写的麻烦。3一 7. 4 算例例土 9设深梁承受巧布载荷，如图3. 36(a)所示,暇定E=l,泊松顿=饥17.不汁 容斗.厚度=】萨为平面应力间萩.|明称取半边结构计沉托构支承.单元划分.节点 编号如图3. 3时S所不.试哽出工一。及)=8 m截面的竖向位移图口 = 3 m截面的 应力分布图，源程序为：!Page149COMMON/X1/NJ,NE,NZ,

2、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 10EO=EO/(1.0-UN*UN)UN=UN/(1.0-UN)10 CALL TOTSTICALL LOADCALL SUPPORCALLSOLVEQCALL STRESSPAUSE!

3、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

4、=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)

5、,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,3DO 10 J=1,6B(I,J)=0.010 CONTINUEB(1,1)=-BJ-BMB(1,3)=BJB(1,5)=BMB(2,2)=-CJ-CMB(2,4)

6、=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,3DO 20 J=1,6B(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,3DO

7、 30 J=1,6S(I,J)=0.0DO 30 K=1,3S(I,J)=S(I,J)+D(I,K)*B(K,J)30 CONTINUEIF(IASK.LE.2) GOTO 50DO 40 I=1,6DO 40 J=1,6EKE(I,J)=0.0DO 40 K=1,3And!*ExchangeS*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(10

8、0,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,NECALL ELEST(MEO,3)DO 30 I=1,3DO 30 II=1,2LH=2*(I-1)+IILDH=2*(JM(MEO,I)-1)+IIDO 30 J=1,3DO 30 JJ=1,2L=2*(J-1)+JJLZ=2*(JM(MEO,J)-1

9、)+JJLD=LZ-LDH+1IF(LD.LE.0) GOTO 30TKZ(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,NJ2P(I)=0.010 CONTINUEIF(NPJ.EQ.0) GOTO

10、30DO 20 I=1,NPJI1=I+1J=IFIX(PJ(I1,2)P(J)=PJ(I1,1)20 CONTINUE30 IF(GAMA.LE.0.0) GOTO 50!Page153DO 40 MEO=1,NECALL ELEST(MEO,1)PE=-GAMA*AE*TE/3.0IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)P(2*IE)=P(2*IE)+PEP(2*JE)=P(2*JE)+PEP(2*ME)=P(2*ME)+PE40 CONTINUERETURNENDSUBROUTINE SUPPORCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND

11、,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,NZMZ=NZC(I)TKZ(MZ,1)=1.0DO 10 J=2,NDDTKZ(MZ,J)=0.010 CONTINUEIF(MZ-NDD)20,20,3020 JO=MZGOTO 4030 JO=NDD40 DO 50 J = 2,JOJ1=MZ-JTKZ(J1+1,J)=0.050 CONTINUEP(MZ)=0.060 CONTINU

12、ERETURNEND!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,NJ1IF(NJ2-K-NDD+1)10,10,2010 IM=NJ2GOTO 3020 IM=K+NDD-130 K1=K+1L=I-K+1C=TKZ(K,L)/TKZ(K,1)LD1=NDD-

13、L+1DO 40 J=1,LD1M=J+I-KTKZ(I,J)=TKZ(I,J)-C*TKZ(K,M)40 CONTINUEP(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=NDDGOTO 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!Page155WR

14、ITE(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 ME

15、O=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

16、) 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,CETA50FORMAT(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) !50FORMAT(4X,2HE=,I3/2X,3HSX=,Fll.3,3X,3HSY=,F11.3,3X

17、,4HTAU=,F11.3/2X,3HSl=,Fll.3,3X,3HS2=,F11.3,3X,4HCET=,F11.3)60 CONTINUERETURNEND输入文件为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,1

18、919,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,1 3,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 TOC o 1-5 h z 0.006.001.006.002.006.003.006.000.005.001.005.002.005.

19、003.005.000.004.001.004.002.004.003.004.000.003.001.003.002.003.003.003.000.002.001.002.002.002.003.002.000.001.001.001.002.001.003.001.000.000.001.000.002.000.003.000.00JD= U=V=1-29766.873-1173917.7502-14003.185-1174018.8753-3753.270-1179518.125407505-26382.471-1072681.5006-10746.993-1

20、073615.0007-2064.593-1082360.750802509-13536.995-964010.125103372.794-970055.125117268.415-989269.125120.000-998401.812137816.581-835383.4381427176.234-861713.9381522063.230-905726.125160.000-927165.1881729514.479-665602.8751853419.637-747340.4381934876.832-839806.812200.000-881219.1252

21、129580.273-416288.7192252944.918-632601.1252317504.195-803765.688240.000-859481.938250.0000.00026-120102.820-583505.37527-76202.375-787347.188280.000-829170.812E= 1SX=-1489.530SY=-101489.383TAU=-1489.531S1=-1467.348S2=-101511.562CET=179.147E=2SX=-1475.844SY=-100654.875TAU=-1790.500S1=-1443.531S2=-10

22、0687.188CET=178.966E=3SX=-7021.670SY=-101597.672TAU=-3741.688S1=-6873.875S2=-101745.469CET=177.738E=4SX=-8067.500SY= -98528.750TAU=-4459.156S1=-7848.227S2= -98748.023CET=177.185E=5SX= -13143.328SY= -99391.750TAU=-1662.500S1=-13111.293S2= -99423.781CET=178.896E=6SX= -14652.781SY= -98337.500TAU=-1501.

23、062S1= -14625.867E= 7SX= -2923.122S1= -2597.762E= 8SX= -716.078S1=0.148E= 9SX= -9188.316S1= -8203.125E= 10SX= -12285.000S1= -10526.887E= 11SX= -14170.516S1= -13801.664E= 12SX= -22797.406S1= -22575.914E= 13SX= -5104.269S1= -4011.727E= 14SX= 969.672S1=5024.582E= 15SX= -14954.572S1= -11552.273E= 16SX=

24、-19890.141S1= -11514.844E= 17SX= -22109.688S1= -20543.453E= 18SX= -35190.453S1= -33280.023E= 19SX= -9785.850S1= -7220.594E= 20SX= 4594.438S1= 20477.398E= 21S2= -98364.414SY=-109168.297 S2=-109493.656SY=-103681.562S2=-104397.789SY=-105121.867S2=-106107.062SY= -95180.250S2= -96938.359SY= -95500.750S2=

25、 -95869.602SY= -91347.000S2= -91568.492SY=-129494.438S2=-130586.977SY=-108176.375S2=-112231.281SY=-110883.469S2=-114285.766SY= -86924.312S2= -95299.609SY= -87301.625S2= -88867.859SY= -77219.000S2= -79129.430SY=-171444.172S2=-174009.422SY=-113592.375S2=-129475.336CET= 178.973TAU= -5888.469CET= 176.83

26、7TAU= -8617.406CET= 175.249TAU= -9771.594CET= 174.243TAU= -12199.594CET= 171.799TAU= -5489.531CET= 176.156TAU= -3902.844CET= 176.752TAU= -11708.750CET= 174.669TAU= -21424.750CET= 169.283TAU= -18383.531CET= 169.515TAU= -25131.188CET= 161.569TAU= -10225.406CET= 171.292TAU= -9162.000CET= 168.222TAU= -2

27、0524.969CET= 172.876SX= -25287.307S1= -16467.512E= 22SX= -30634.422S1= -1543.715E= 23SX= -34259.609S1= -29220.699E= 24SX= -43958.047S1= -30369.627E= 25SX= -19028.160S1= -13907.055E= 26SX= 3973.812S1= 54459.047E= 27SX= -39180.809S1= -23409.074E= 28SX= -42804.766S1= 22662.211E= 29SX= -42224.094S1= -32436.225E= 30SX= -21830.422S1=239.594E= 31SX= -48815.199S1= -32570.844E= 32SX=-132271.750S1= 76087.781E= 33SX= -45090.102S1= -45034.867E= 34SX= 42332.711S1= 70218.328E= 35SX= -208