matlab实体单元,基于matlab三维实体六单元的分析

基于matlab三维实体六单元的分析

Matlab 三维实体六单元的分析摘要：Matlb 三维实体六单元 分析目标：通过对三维实体的分析，进一步掌握有限元计算的方法和步骤，了解实体上各点的应变状态。模型：如图 1 所示的一个长方体，长宽高分别为 40cm，10cm，30齿面，分成六个单元，每个单元有 8 个节点，且几何尺寸相同，都为 40cm X 10cm X 5cm。长方体顶部及底部为固定约束，在节点 10处沿 z 方向作用 F=1000N，E=2.06 X 1011N/m， 。 (老师图我不会画0.3就没画了)内容：Matlab 程序如下：clearclc% 3D-solid elments for a shaft%% (x1 y1 z1)at shortE=2.06e11; Nu=0.3;h=0.04;%(1)球形函数矩阵x1=0;y1=0;z1=0;x2=0;y2=0;z2=0.04;x3=0.01;y3=0;z3=0.04;x4=0.01;y4=0;z4=0;x5=0;y5=0.005;z5=0;x6=0;y6=0.005;z6=0.04;x7=0.01;y7=0.005;z7=0.04;x8=0.01;y8=0.005;z8=0;x=sym( x );y=sym( y );z=sym( z );v=[1 x y z x*y y*z z*x x*y*z];m=.[1 x1 y1 z1 x1*y1 y1*z1 z1*x1 x1*y1*z11 x2 y2 z2 x2*y2 y2*z2 z2*x2 x2*y2*z21 x3 y3 z3 x3*y3 y3*z3 z3*x3 x3*y3*z31 x4 y4 z4 x4*y4 y4*z4 z4*x4 x4*y4*z41 x5 y5 z5 x5*y5 y5*z5 z5*x5 x5*y5*z51 x6 y6 z6 x6*y6 y6*z6 z6*x6 x6*y6*z61 x7 y7 z7 x7*y7 y7*z7 z7*x7 x7*y7*z71 x8 y8 z8 x8*y8 y8*z8 z8*x8 x8*y8*z8];mm=inv(m);N=v*mm;%u={u1 v1 w1 u2 v2 w2 u3 v3 w3,u8 v8 w8 } uN=[N(1) 0 0 N(2) 0 0 N(3) 0 0 N(4) 0 0 N(5) 0 0 N(6) 0 0 N(7) 0 0 N(8) 0 0];vN=[0 N(1) 0 0 N(2) 0 0 N(3) 0 0 N(4) 0 0 N(5) 0 0 N(6) 0 0 N(7) 0 0 N(8) 0];wN=[0 0 N(1) 0 0 N(2) 0 0 N(3) 0 0 N(4) 0 0 N(5) 0 0 N(6) 0 0 N(7) 0 0 N(8)];% 2 求单元的刚度矩阵B1=diff(uN,x);B2=diff(vN,y);B3=diff(wN,z);B4=diff(uN,y)+diff(vN,x);B5=diff(vN,z)+diff(wN,y);B6=diff(wN,x)+diff(uN,z);B=[B1;B2;B3;B4;B5;B6];BT=transpose(B);D=.E/[(1+Nu)*(1-2*Nu)]*[1-Nu Nu Nu 0 0 0;Nu 1-Nu Nu 0 0 0;Nu Nu 1-Nu 0 0 0;0 0 0 (0.5-Nu) 0 0;0 0 0 0 (0.5-Nu) 0;0 0 0 0 0 (0.5-Nu)];K=BT*D*B;KI=int(K,x,x1,x7);KI=int(KI,y,y1,y7);KI=int(KI,z,z1,z7);KI=vpa(KI);Ke1=KI;Ke2=KI;Ke3=KI;Ke4=KI;Ke5=KI;Ke6=KI;% 3 应用转化矩阵法进行整体刚度矩阵的组集Ndof=3*4*7;G1=zeros(24,Ndof);G1(:,12*0+1:12*0+24)=eye(24,24);G2=zeros(24,Ndof);G2(:,12*1+1:12*1+24)=eye(24,24);G3=zeros(24,Ndof);G3(:,12*2+1:12*2+24)=eye(24,24);G4=zeros(24,Ndof);G4(:,12*3+1:12*3+24)=eye(24,24);G5=zeros(24,Ndof);G5(:,12*4+1:12*4+24)=eye(24,24);G6=zeros(24,Ndof);G6(:,12*5+1:12*5+24)=eye(24,24);K1=G1 *Ke1*G1;K2=G2 *Ke2*G2;K3=G3 *Ke3*G3;K4=G4 *Ke4*G4;K5=G5 *Ke5*G5;K6=G6 *Ke6*G6;K=K1+K2+K3+K4+K5+K6;%4 引入约束条件F=zeros(Ndof,1);F(10*3)=1000;KX=K(13:Ndof-12,13:Ndof-12);FX=F(13:Ndof-12,1);% 5 进行求解u=inv(KX)*FX;u=double(u)求得各节点的位移如下u =1.0e-005 *5: -0.00640.0001-0.00026: 0.01020.00240.02907: 0.0078-0.00430.01628: -0.00910.0031-0.00159: -0.0124-0.0019-0.025810: 0.01980.00100.100211： 0.0136-0.00210.003812：-0.01590.00130.008813： -0.0130-0.00350.004014： 0.0178-0.00130.038515： 0.01370.00220.023016：-0.0148-0.00260.000417： -0.0091-0.00180.007018： 0.0113-0.00200.019919：0.00950.00440.014520：-0.0108-0.00550.002121：-0.0034-0.00030.004422： 0.0044-0.00160.008823： 0.00450.00370.006324： -0.0057- 0.00460.0017结果分析：(1) 把上述计算数据和书上划分三单元的计算结果比较，数据相差很大，看来三维实体的计算精度和单元的划分有很大关系。(2) 10 节点处的 Z 方向的位移最大 ,10 节点附近5max0.12wm的 11、14 节点各方向的位移也较其它节点大。参考文献：[1]韩清凯，孙伟。弹性力学及有限元法基础教程。沈阳：东北大学出版社，2009 。