湖南科技大学有限元课程设计

1、湖南科技大学程设计报告有限元法理论与应用课程设计平面应力问题有限元分析土木工程学院专业班级2013级工程力学2班1302060213谢献忠学生姓名 学生学号 指导教师零一六年六月. 课程设计题目.有限元计算模型三. 基本理论四. 计算结果及分析五. 源程序目录有限元法理论与应用课程设计正文.课程设计题目 我选择的题目是：题3 图示为一带圆孔的单位厚度(1M的正方形平板，在X方向作用均布压力0.25MPa试分别用不同数量的三节点常应变单元对平板进行有限元分析。在丫轴上，孔边应力的精确解为：6 x=-0.75MPa,在X轴上，孔边应力的精确解为：6 y=0.25MPa0. 25MPa0. 25MP

2、a6mm48m分为180个单元,112个节点。.有限元计算模型 (1)该平板的受力属于平面应力问题，模型受到对称外荷载的作用，其实际受力情况也是 对称的。(2)对该平板可以简化为一平板倒圆角的模型(即取平板的四分之一部分) ，在对称荷载作 用下平板的内部受力图也应该是对称的。且平板只受到沿 X方向的外荷载作用。板面上无 作用力，故Z方向没有正应力作用。(3)外荷载的处理是将均布荷载简化为作用在右边界处的节点承受集中荷载，并通过叠加 作用分布在各个边界节点处。(4)简化后的平板可以看做倒圆左边界视作X约束而丫无约束，倒圆下边界方向视有丫约束X无约束，并且用三节点划分单元网格。 因为中间圆孔存在应

3、力集中，所以在孔的附近需要将网格细化，远离网格的区域则不 需要画细。网格的划分图(因节点数太多在孔边处集中难以放大，此处划分的图用简图 )9三.基本理论单元刚度矩阵Ke = ff BTDBtcXdy，总体刚度矩阵K=刀L(i)Tk(i)L(i)弹性矩阵D= E(1-u(1+u)(1-2u) 1u/(1-u)u/(1-u)孔边单元的应力u/(1-u)坐标变换 u= a 1+ a 2E + a 3 n + a 4E n ,v=a 5+ a 6E + a 7 n + a 8E n荷载等效移置Re=NTcG+ 职Tqtds+yiN|T ptdxdy四.计算结果及分析孔边单元的应力123132eiSE1

4、-o.osmeC.001932-0. 10T35-0. 000064-0.E29811-0. 15395C. 166320.0634550.0969250.0550310.0036424041968-0 01BO30 0S03510.00912320.0721860.070470 067399111眈11L22151152-0 3eT9S-0 444-0.52571-0. 3255657435-0.39497-0.067543-0.027196-0. 0953閃-0.062294-0.061221-0.05990 11590 100530.10S2&0. i04se0 04312T0.0912

5、37孔边单元的应变1313EG16-2.59416-005-5.7015&-009-4.44TTe-aoe-3 2291e-0a8-T.695Se-00e-5.4954-OOe5 8304-008. 0958e-00e3 9TU-OOS2.6633-0082.4195e-00e2. 24554-006-6.4344e-0092 154&-ooe3.S5Te-0092 57Sle-OOS2.514Se-00e3.1354e-00S9102 11IZlIZZ15115-1.257e-007-7.e9Se-ooe-1.657t-00T-1.0229 #-007-l.e533e-007-1.2567e

6、-007l.6281e-00e1 5429e-ooeZ.OTS3e-OOe1. lT91e-00e3. 70284-0081.953e-00e4. 1361e-QOS3.53T4e-00e3. 3633&-0083.749e-00e1.539e-OOS3.5e&-ooe位移图:孔边单元的应力12 129305753esee1-25日4日e-5.7034 e-. J-4 449e-008-32302&-7.699e-008-5.497e-008-1.22Be-0D7-7 895e-DDB5.S412e-0082.0973eOTi3 9731e-0Dg2.6699e-0DB2.4209 e-D08

7、2.246Be-ra:1.629e-0Q81.5438e-008-6.4385e-.2.1548eOTi3.2519e-0092.5767e-OOB2 5l5e-0083.1371 e-0084.13686-0083.58B2e-OOB115116142143-1.2345e-.,.-9 34390-.,-1.257e-D0713 241 一2.0373e-ra2.2715eOT1.95370-0081.4164e-Q0e-I.OOtSe-.17.9497e-Q091 13.2589 e-008-I.OIBe-OOB位移图（180单元105节点）孔边单元的应变孔边单元的应力1 123334&5

8、EG-0 02777610.0019356-0.10736-0 0B0D73-0.22984-0.158980.166360.063476QDB69510.056D490 0036532D. 019691-0 0180360,060370,00911840 0721940,0704730.0878989697I1291130159160-025-.525770.3256-0.3253-0 57527-3.1282e-.-0.0E7551I-0.095374-0 06229867776-0.362534 0S19e-0090.11591.10B270.10490 0352330 54541孔边单

9、元的应变123334esGS1-2.5944 e-OOB-5.7024e-009-4.44B3e-.-3.2296e-.-7E978e-D0B-5.49G1e-0D85S397e-D0S2.D965e-0QS3 972100日2.669e-0082.42D2e-D0B2.2461 e-DDB6 43E3e-0092,1544e-00B3.254e-D092.57B3e-ra2.5l49e-OOS3 1367e-00896日丁I129130159I160-8.33330-008-1.2259e-007 1.6572 e-007-1.D23e-DD7-4.0654e-008M .555e-D072

10、.5e-0Q81. 2859-0082.0786e-0081.1794G-OOB-1.9340-007-6 3317e-,1.45G7e-D154.1364e-0D83 86373.7433e1.2573e-0Da1.94G4e-0D7位移图结论：对比不同节点数量的方案，可得不同的方案下求得的孔边应力应变大小接近，且单 元数多的方案，其应力应变略大于单元数少的方案。由于变形较小，在同放大500000倍的情况下位移图基本一致，同时单元数越多越接近精确解但小于精确解。五源程序clear clc n=0;RR=3 4 6 8 10 15 25 50 100 300 400 2000 5000 100

11、00 15000 20000 0; for i=0:15:90for j=1:16 n=n+1;co1(n,:)=RR(j)*cos(i*pi/180) RR(j)*sin(i*pi/180);endendx=24,0;24,24*tan(15*pi/180);24,24*tan(30*pi/180);24,24*tan(45*pi/180);24*tan(30*pi/180),24;24*tan(15*pi/180),24;0,24*1000;for j=1:7co1(j*16,:)=x(j,:);end% 生成节点坐标矩阵clear i j n x RR% n=2;n1=1;n2=2;n3

12、=18;for i=1:6for j=1:15 ijm(n,:)=n1 n2 n3;n1=n1+1; n2=n2+1;n3=n3+1;n=n+2;endn1=n1+1; n2=n2+1;n3=n3+1;end%生成偶数单元的ijm矩阵clear n n1 n2 n3 j i n=1;n1=1;n2=18;n3=17;for i=1:6for j=1:15 ijm(n,:)=n1 n2 n3;n1=n1+1; n2=n2+1;n3=n3+1;n=n+2;endn1=n1+1; n2=n2+1;n3=n3+1;end%生成奇数单元的ijm矩阵% noc=180;%单元数nod=112%节点数 do

13、f=2*nod;% 自由度数 EE=3e6;%弹性模量vv=0.3;% 泊松比 %for i=1:6;for j=1:30;E(i-1)*30+j)=EE;v(i-1)*30+j)=vv;endendclear i j EE vv%插入各单元弹性模量、泊松比%l=48000;% 长度h=48000% 高度 t=1000;% 厚度 %(二) 计算单元刚度矩阵并组集整体刚度矩阵K=zeros(2*nod,2*nod);% 整体刚度矩阵赋零初值for i=1:nocmA=1,co1(ijm(i,1),1) co1(ijm(i,1),2); 1,co1(ijm(i,2),1) co1(ijm(i,2)

14、,2); 1,co1(ijm(i,3),1) co1(ijm(i,3),2);% aa(i)=0.5*det(mA);A(i)=aa(i);% 单元面积clear aa单元节点坐标矩阵ma=co1(ijm(i,2),1) co1(ijm(i,2),2); co1(ijm(i,3),1) co1(ijm(i,3),2); a(i,1)=det(ma);ma=co1(ijm(i,3),1) co1(ijm(i,3),2); co1(ijm(i,1),1) co1(ijm(i,1),2);a(i,2)=det(ma);ma=co1(ijm(i,1),1) co1(ijm(i,1),2); co1(i

15、jm(i,2),1) co1(ijm(i,2),2);a(i,3)=det(ma);%生成参数 a 矩阵mb=1,co1(ijm(i,2),2);1,co1(ijm(i,3),2);b(i,1)=-det(mb);mb=1,co1(ijm(i,3),2);1,co1(ijm(i,1),2);b(i,2)=-det(mb);mb=1,co1(ijm(i,1),2);1,co1(ijm(i,2),2);b(i,3)=-det(mb);%生成参数 b 矩阵mc=1,co1(ijm(i,2),1);1,co1(ijm(i,3),1);c(i,1)=det(mc);mc=1,co1(ijm(i,3),1

16、);1,co1(ijm(i,1),1);c(i,2)=det(mc);mc=1,co1(ijm(i,1),1);生成参数 c 矩阵1,co1(ijm(i,2),1);c(i,3)=det(mc);%clear mb mc mA ma end for i=1:nocmB=b(i,1),0;0,c(i,1); c(i,1),b(i,1); mB=0.5/A(i)*mB;B(:,1:2,i)=mB;mB=b(i,2),0;0,c(i,2); c(i,2),b(i,2); mB=0.5/A(i)*mB;B(:,3:4,i)=mB;mB=b(i,3),0;0,c(i,3); c(i,3),b(i,3);

17、 mB=0.5/A(i)*mB;%11B(:,5:6,i)=mB;%生成应变矩阵 B（3,6, 单元号 ）mD=1,v(i),0;v(i),1,0;0,0,(1-v(i)/2);D(:,:,i)=E(i)/(1-v(i)*v(i)*mD;% end平面应力问题弹性矩阵 D（3,3 ，单元号）clear i mB mDfor i=1:nocKe(:,:,i)=B(:,:,i)*D(:,:,i)*B(:,:,i)*t*A(i);%L=zeros(6,2*nod);L(1:2,(ijm(i,1)-1)*2+1:ijm(i,1)*2)=eye(2,2);%1L(3:4,(ijm(i,2)-1)*2+1

18、:ijm(i,2)*2)=eye(2,2);%2L(5:6,(ijm(i,3)-1)*2+1:ijm(i,3)*2)=eye(2,2);%K=K+L*Ke(:,:,i)*L;%end单元刚度矩阵 k（6,6, 单元号 ）生成定位矩阵 L 将单刚矩阵 k 组集成总刚矩阵 Kclear i0%(三) 形成节点荷载 R=zeros(2*nod,1); for i=1:3mR=zeros(2*nod,1);mR(16*i*2-1,1)=1;mR(2*(i+1)*16-1,1)=1;xR=zeros(2*nod,1);ll=co1(i+1)*16,2)-co1(i*16,2);xR=-0.25*1000

19、*ll/2.*mR;R=xR+R;end clear xR mR iR0=112 111 110 109 108 107 106 105 104 103 102 101 100 99 98 97;R1=16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1;for i=1:16R(R0(i)*2-1,:)=;K(R0(i)*2-1,:)=;K(:,R0(i)*2-1)=;end clear i for i=1:16R(R1(i)*2,:)=;K(R1(i)*2,:)=;K(:,R1(i)*2)=;end%上面带都是删除位移在X方向或丫方向上位移为0的节点所对应的总纲矩阵上

20、的行和列% uuww=inv(K)*R;% 计算位移 z=zeros(2,1,112);n=1;for i=1:16z(:,:,i)=uuww(n,1);0;n=n+1;end clear i for i=17:96z(:,:,i)=uuww(n,1);uuww(n+1,1);n=n+2;end clear i for i=97:112z(:,:,i)=0;uuww(n,1);n=n+1;end for i=1:112uw(2*i-1:2*i,1)=z(:,:,i);end%还原删除位移在 X方向或Y方向上位移为0的节点所对应的总纲矩阵上的行和列 clear uuww z i整体节点水平位移向

21、量整体节点垂直位移向量u=zeros(nod,1);w=zeros(nod,1); u(1:nod,1)=uw(1:2:2*nod-1);% w(1:nod,1)=uw(2:2:2*nod);% det=zeros(6,noc);for i=1:nocdet(1,i)=u(ijm(i,1);det(2,i)=w(ijm(i,1);det(3,i)=u(ijm(i,2);det(4,i)=w(ijm(i,2);det(5,i)=u(ijm(i,3);det(6,i)=w(ijm(i,3);endfor i=1:nocepsl(:,i)=B(:,:,i)*det(:,i); %计算每一个单元的应变

22、14end计算每一个单元的应力for i=1:nocsigm(:,i)=D(:,:,i)*epsl(:,i); %end uu=u*500000;ww=w*500000;% 水平、垂直位移放大 500,000 倍co2(:,1)=co1(:,1)+uu;co2(:,2)=co1(:,2)+ww;% 变形后的 X、Y 坐标 y=co1(1,:);co1(2,:);co1(3,:);co1(4,:);co1(5,:);co1(6,:);co1(7,:);co1(8,:);co1(9,:);co1(10,:) ;co1(11,:);co1(12,:);co1(13,:);co1(14,:);co1(

23、15,:);co1(16,:);co1(32,:);co1(48,:);co1(64,:);co 1(80,:);co1(96,:);co1(112,:);co1(111,:);co1(110,:);co1(109,:);co1(108,:);co1(107,:);co1(106,:7,:);co1(81,:);co1(65,:);co1(49,:);co1(33,:);co1(17,:);co1(1,:); % yy=co2(1,:);co2(2,:);co2(3,:);co2(4,:);co2(5,:);co2(6,:);co2(7,:);co2(8,:);co2(9,:);co2(10,: ););co1(105,:);co1(104,:);co1(103,:);co1(102,:);co1(101,:);co1(100,:);co1(99,:);co1(98,:);co1(9变形前的轮廓线co2(11,:);co2