下载文档原格式
计算力学课程设计报告有限元程序求解弹性力学平面问题专业：班级：姓名：学号：指导教师：有限元程序求解弹性力学平面问题设计目的:1、学习有限元程序求解弹性力学平面问题的方法；2、学习有限元程序编写技巧；3、加深对有限元方法的理解；4、锻炼处理复杂弹性力学问题的能力。
题一：例3.9设深梁承受均布荷载，如下图（a）所示。
假定E=1，泊松比=0.17μ，不计容重，厚度t=1m，为平面应力问题。
因对称去半边结构进行计算，结构支承、单元划分、节点编号如图（b）所示。
试画出y=0及y=6m截面的竖向位移图，x=3m截面的xσ应力分布图。
50kN100kN100kN50kNy1、有限元Fortran源程序如下：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)CALL DATAIF(IND.EQ.0)GOTO 10EO=EO/(1.0-UN*UN)UN=UN/(1.0-UN)10 CALL TOTSTICALL LOADCALL SUPPORCALL SOLVEQCALL STRESSSTOPENDSUBROUTINE DATACOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/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)PRINT*,"INPUT:NJ,NE,NZ,NDD,NPJ,IND"READ(5,*)NJ,NE,NZ,NDD,NPJ,INDNJ2=NJ*2NPJ1=NPJ+1PRINT*,"INPUT:EO,UN,GAMA,TE"READ(5,*)EO,UN,GAMA,TEPRINT*,"INPUT:JM"READ(5,*)((JM(I,J),J=1,3),I=1,NE)PRINT*,"INPUT:CJZ"READ(5,*)((CJZ(I,J),J=1,2),I=1,NJ)PRINT*,"INPUT:NZC"READ(5,*)(NZC(I),I=1,NZ)PRINT*,"INPUT:PJ"READ(5,*)((PJ(I,J),J=1,2),I=1,NPJ1)OPEN(100,FILE='1.TXT')WRITE(100,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,AE COMMON/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,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)=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)B(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 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,6DO 40 K=1,3EKE(I,J)=EKE(I,J)+S(K,I)*B(K,J)*AE*TE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE TOTSTICOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/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 20 I=1,NJ2DO 20 J=1,NDDTKZ(I,J)=0.020 CONTINUEDO 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)+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,AE COMMON/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 30DO 20 I=1,NPJI1=I+1J=IFIX(PJ(I1,2))20 CONTINUE30 IF(GAMA.LE.0.0) GOTO 50DO 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 CONTINUE50 CONTINUERETURNENDSUBROUTINE SUPPORCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/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 CONTINUERETURNENDSUBROUTINE SOLVEQCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/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 50 K=1,NJ1IF(NJ2-K-NDD+1)10,10,2010 IM=NJ2GOTO 3020 IM=K+NDD-130 K1=K+1DO 50 I=K1,IML=I-K+1C=TKZ(K,L)/TKZ(K,1)LD1=NDD-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-I1IF(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 CONTINUEOPEN(101,FILE='2.TXT')WRITE(101,110)(I,P(2*I-1),P(2*I),I=1,NJ)110 FORMAT(2X,3HJD=,5X,2HU=,12X,2HV=/(I4,3X,F11.3,2X,F11.3)) RETURNENDSUBROUTINE STRESSCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/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,3LH=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 4030 CETA=0.0OPEN(102,FILE='3.TXT')40 WRITE(102,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)60 CONTINUERETURNEND2、输入数据：28,36,9,10,4,0 1,0.17,0,1 1,5,22,5,62,6,33,6,73,7,44,7,85,9,66,9,106,10,7 7,10,11 7,11,8 8,11,12 9,13,10 10,13,14 10,14,11 11,14,15 11,15,12 12,15,16 13,17,14 14,17,1814,18,15 15,18,19 15,19,16 16,19,20 17,21,18 18,21,22 18,22,19 19,22,23 19,23,20 20,23,24 21,25,22 22,25,26 22,26,23 23,26,27 23,27,24 24,27,28 0,61,62,63,60,51,52,53,50,4 1,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,55 0,0-50000,2-100000,4-100000,6-50000,83、输出结果：⑴、节点坐标：NO X Y1 0.00 6.002 1.00 6.003 2.00 6.004 3.00 6.005 0.00 5.006 1.00 5.007 2.00 5.008 3.00 5.009 0.00 4.0010 1.00 4.0011 2.00 4.0012 3.00 4.0013 0.00 3.0014 1.00 3.0015 2.00 3.0016 3.00 3.0017 0.00 2.0018 1.00 2.0019 2.00 2.0020 3.00 2.0021 0.00 1.0022 1.00 1.0023 2.00 1.0024 3.00 1.0025 0.00 0.0026 1.00 0.0027 2.00 0.0028 3.00 0.00⑵、节点位移（m）：JD= U= V=1 -29767.031 -1173919.8752 -14003.248 -1174020.7503 -3753.289 -1179520.1254 0.000 -1181721.8755 -26382.504 -1072683.5006 -10746.987 -1073616.8757 -2064.578 -1082362.6258 0.000 -1085875.3759 -13536.982 -964012.00010 3372.830 -970057.00011 7268.439 -989270.93812 0.000 -998403.81213 7816.650 -835385.12514 27176.352 -861715.688 15 22063.291 -905727.93816 0.000 -927167.25017 29514.604 -665604.25018 53419.770 -747342.18819 34876.914 -839808.62520 0.000 -881221.00021 29580.316 -416289.65622 52945.020 -632602.68823 17504.215 -803767.50024 0.000 -859483.75025 0.000 0.00026 -120103.062 -583506.75027 -76202.484 -787348.93828 0.000 -829172.562⑶、单元应力（Pa）：E= 1SX= -1489.449 SY=-101489.578 TAU= -1489.493 S1= -1467.270 S2=-101511.766 CET= 179.147 E= 2SX= -1475.786 SY=-100654.781 TAU= -1790.448 S1= -1443.477 S2=-100687.094 CET= 178.966 E= 3SX= -7021.617 SY=-101597.578 TAU= -3741.747 S1= -6873.812 S2=-101745.383 CET= 177.738 E= 4SX= -8067.512 SY= -98528.883 TAU= -4459.167 S1= -7848.234 S2= -98748.156 CET= 177.185E= 5SX= -13143.324 SY= -99391.883 TAU= -1662.572 S1= -13111.285 S2= -99423.922 CET= 178.896E= 6SX= -14652.791 SY= -98337.594 TAU= -1501.185 S1= -14625.871 S2= -98364.516 CET= 178.973E= 7SX= -2923.120 SY=-109168.492 TAU= -5888.416 S1= -2597.762 S2=-109493.844 CET= 176.837 E= 8SX= -716.068 SY=-103681.609 TAU= -8617.434 S1= 0.160 S2=-104397.844 CET= 175.249 E= 9SX= -9188.308 SY=-105121.906 TAU= -9771.604 S1= -8203.109 S2=-106107.102 CET= 174.243E= 10SX= -12285.007 SY= -95180.242 TAU= -12199.550S1= -10526.906 S2= -96938.344 CET= 171.799 E= 11SX= -14170.538 SY= -95500.742 TAU= -5489.623 S1= -13801.672 S2= -95869.609 CET= 176.156 E= 12SX= -22797.447 SY= -91347.070 TAU= -3902.925 S1= -22575.945 S2= -91568.570 CET= 176.752 E= 13SX= -5104.274 SY=-129494.578 TAU= -11708.809 S1= -4011.730 S2=-130587.125 CET= 174.669 E= 14SX= 969.693 SY=-108176.422 TAU= -21424.820 S1= 5024.629 S2=-112231.359 CET= 169.283 E= 15SX= -14954.604 SY=-110883.578 TAU= -18383.518 S1= -11552.312 S2=-114285.867 CET= 169.515 E= 16SX= -19890.213 SY= -86924.367 TAU= -25131.225S1= -11514.891 S2= -95299.688 CET= 161.569 E= 17SX= -22109.713 SY= -87301.680 TAU= -10225.519S1= -20543.445 S2= -88867.945 CET= 171.292 E= 18SX= -35190.500 SY= -77218.953 TAU= -9162.105 S1= -33280.023 S2= -79129.430 CET= 168.221 E= 19SX= -9785.861 SY=-171444.500 TAU= -20525.008 S1= -7220.609 S2=-174009.750 CET= 172.876 E= 20SX= 4594.449 SY=-113592.445 TAU= -46145.883 S1= 20477.516 S2=-129475.516 CET= 161.007 E= 21SX= -25287.357 SY=-118672.375 TAU= -30023.787 S1= -16467.543 S2=-127492.188 CET= 163.629 E= 22SX= -30634.465 SY= -71127.195 TAU= -44991.484 S1= -1543.734 S2=-100217.922 CET= 147.114 E= 23SX= -34259.684 SY= -71743.445 TAU= -14638.012 S1= -29220.699 S2= -76782.422 CET= 161.005 E= 24SX= -43958.176 SY= -53419.129 TAU= -17697.580S1= -30369.762 S2= -67007.547 CET= 142.483 E= 25SX= -19028.232 SY=-252549.391 TAU= -34958.820 S1= -13907.102 S2=-257670.531 CET= 171.666 E= 26SX= 3973.836 SY=-114063.977 TAU= -92238.578 S1= 54459.203 S2=-164549.344 CET= 151.307 E= 27SX= -39180.902 SY=-121400.266 TAU= -39312.672 S1= -23409.199 S2=-137171.969 CET= 158.140 E= 28SX= -42804.859 SY= -43317.969 TAU= -65723.141S1= 22662.227 S2=-108785.055 CET= 135.112 E= 29SX= -42224.188 SY= -43219.219 TAU= -10273.362S1= -32436.301 S2= -53007.105 CET= 136.386 E= 30SX= -21830.449 SY= -25448.408 TAU= -23810.377S1= 239.566 S2= -47518.426 CET= 137.172 E= 31SX= -48815.301 SY=-424588.281 TAU= -79800.312 S1= -32570.906 S2=-440832.688 CET= 168.494 E= 32SX=-132272.047 SY= -71582.141 TAU=-175409.688 S1= 76088.000 S2=-279942.188 CET= 130.093 E= 33SX= -45090.219 SY= -56761.309 TAU= 804.799 S1= -45034.984 S2= -56816.547 CET= 3.926 E= 34SX= 42332.844 SY= -9222.024 TAU= -47066.441 S1= 70218.484 S2= -37107.668 CET= 149.354 E= 35SX= -20899.361 SY= -19971.461 TAU= 16235.217 S1= -4193.565 S2= -36677.254 CET= 45.818 E= 36SX= 73164.031 SY= -17873.285 TAU= -17873.346 S1= 76547.367 S2= -21256.619 CET= 169.2814、数据处理：由以上运行结果，并结合对称特点可得y=0及 y=6面上节点竖向位移值如下表：根据上表位移值，可作出y=0及 y=6截面的竖向位移图如下：整理x=3m 处各点的x σ值时，图（b ）中的b 、c 、d 、e 、f 、g 六点按二单元平均法整理；a 、h 两点为边界点，在a 点附近应力x σ变化显著，用四点插值公式计算a 点的应力；而h 点用三点插值公式计算应力。
2.6 四结点四边形单元(The four-node quadrilateral element)前面介绍了四结点的矩形单元其位移函数:xy y x U 4321αααα+++=xy y x V8765αααα+++=为双线性函数，应力，应变在单元内呈线性变化，比常应力三角形单元精度高。
但它对边界要求严格。
本节介绍的四结点四边形等参元，它不但具有较高的精度，而且其网格划分也不受边界的影响。
对任意四边形单元（图见下面）若仍直接采用前面矩形单元的位移函数，在边界上它便不再是线性的（因边界不与x,y 轴一致），这样会使得相邻两单元在公共边界上的位移可能会出现不连续现象（非协调元），而使收敛性受到影响。
可以验证，利用坐标变换就能解决这个问题，即可以通过坐标变换将整体坐标中的四边形（图a ）变换成在局部坐标系中与四边形方向无关的边长为2的正方形。
正方形四个结点i,j,m,p 按反时钟顺序对应四边形的四个结点i j m p 。
正方形的 1-=η 和 1=η 二条边界，分别对应四边形的i ，j 边界和p,m 边界；ξ=-1和ξ=+1分别对应四边形的i ，p 边界和j ，m 边界。
如果用二组直线等分四边形的四个边界线段，使四边形绘成一个非正交网格，那么该非正交网格在正方形上对应着一个等距离的规则网格（见图a, b ）。
当然, 局部坐标上的A 点与整体坐标的A 点对应。
一、四结点四边形等参单元的形函数及坐标变换由于可以将整体坐标下的四边形单元变换成局部坐标下的正方形单元，对于这种正方形单元，自然仍取形函数为： ξηαηαξαα2321+++=U ξηαηαξαα8765+++=V引入边界条件，即可得位移函数：∑=ijmpi i U N Ui ijmpi V N V ∑==写成矩阵形式：{}{}[]{}ee p i p i ed N d N N N N V U f =⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧=000 式中形函数: ()()()ηηξξηξi i i N ++=1141, ()p m j i ,,, 按照等参元的定义,我们将坐标变换式亦取为: p p m m j j i i i ijmpi x N x N x N x N x N x +++==∑p p m m j j i i i ijmpi y N y N y N y N y N y +++==∑ ()162-- 式中形函数N 与位移函数中的完全一致。
