有限元课程设计

一．问题描述如图所示的平面矩形结构，设E=1，NU=0.25，h=1，考虑以下约束和外载：位移边界条件BC（u）：UA=0，VA=0，UD=0，力边界条件BC（p）：在CD边上有均布载荷q=1，建模情形：使用四个四节点矩形单元，试在该建模情形下，求各节点的位移以及各个单元的应力分布。二．Matlab程序(1).函数定义：functionk=Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)symsst;a=(yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;B1=[a*(t-1)/4-b*(s-1)/40;0c*(s-1)/4-d*(t-1)/4;c*(s-1)/4-d*(t-1)/4a*(t-1)/4-b*(s-1)/4];B2=[a*(1-t)/4-b*(-1-s)/40;0c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4a*(1-t)/4-b*(-1-s)/4];B3=[a*(t+1)/4-b*(s+1)/40;0c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4a*(t+1)/4-b*(s+1)/4];B4=[a*(-1-t)/4-b*(1-s)/40;0c*(1-s)/4-d*(-1-t)/4;c*(1-s)/4-d*(-1-t)/4a*(-1-t)/4-b*(1-s)/4];Bfirst=[B1B2B3B4];Jfirst=[01-tt-ss-1;t-10s+1-s-t;s-t-s-10t+1;1-ss+t-t-10];J=[xixjxmxp]*Jfirst*[yi;yj;ym;yp]/8;B=Bfirst/J;ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endBD=J*transpose(B)*D*B;r=int(int(BD,t,-1,1),s,-1,1);z=h*r;k=double(z);endfunctionz=Quad2D4Node_Assembly(KK,k,i,j,m,p)DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;DOF(5)=2*m-1;DOF(6)=2*m;DOF(7)=2*p-1;DOF(8)=2*p;forn1=1:8forn2=1:8KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;endfunctionstress=Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)symsst;a=(yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;B1=[a*(t-1)/4-b*(s-1)/40;0c*(s-1)/4-d*(t-1)/4;c*(s-1)/4-d*(t-1)/4a*(t-1)/4-b*(s-1)/4];B2=[a*(1-t)/4-b*(-1-s)/40;0c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4a*(1-t)/4-b*(-1-s)/4];B3=[a*(t+1)/4-b*(s+1)/40;0c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4a*(t+1)/4-b*(s+1)/4];B4=[a*(-1-t)/4-b*(1-s)/40;0c*(1-s)/4-d*(-1-t)/4;c*(1-s)/4-d*(-1-t)/4a*(-1-t)/4-b*(1-s)/4];Bfirst=[B1B2B3B4];Jfirst=[01-tt-ss-1;t-10s+1-s-t;s-t-s-10t+1;1-ss+t-t-10];J=[xixjxmxp]*Jfirst*[yi;yj;ym;yp]/8;B=Bfirst/J;ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endstr1=D*B*u;str2=subs(str1,{s,t},{0,0});stress=double(str2);end（2）.计算部分E=1;NU=0.25;h=1;ID=1;k1=Quad2D4Node_Stiffness(E,NU,h,1,1,0.5,1,0.5,0.5,1,0.5,ID);k2=Quad2D4Node_Stiffness(E,NU,h,1,0.5,0.5,0.5,0.5,0,1,0,ID);k3=Quad2D4Node_Stiffness(E,NU,h,0.5,1,0,1,0,0.5,0.5,0.5,ID);k4=Quad2D4Node_Stiffness(E,NU,h,0.5,0.5,0,0.5,0,0,0.5,0,ID);KK=zeros(18,18);KK=Quad2D4Node_Assembly(KK,k1,1,6,5,2);KK=Quad2D4Node_Assembly(KK,k2,2,5,4,3);KK=Quad2D4Node_Assembly(KK,k3,6,7,8,5);KK=Quad2D4Node_Assembly(KK,k4,5,8,9,4)k=KK([1:12,14:16],[1:12,14:16]);p=[0;-0.25;0;0;0;0;0;0;0;0;0;-0.5;-0.25;0;0];u=k\pU=[u(1:12);0;u(13:15);0;0];u1=[U(1);U(2);U(11);U(12);U(9);U(10);U(3);U(4)];stress1=Quad2D4Node_Stress(E,NU,1,1,0.5,1,0.5,0.5,1,0.5,u1,ID)u2=[U(3);U(4);U(9);U(10);U(7);U(8);U(5);U(6)];stress2=Quad2D4Node_Stress(E,NU,1,0.5,0.5,0.5,0.5,0,1,0,u2,ID)u3=[U(11);U(12);U(13);U(14);U(15);U(16);U(9);U(10)];stress3=Quad2D4Node_Stress(E,NU,0.5,1,0,1,0,0.5,0.5,0.5,u3,ID)u4=[U(9);U(10);U(15);U(16);U(17);U(18);U(7);U(8)];stress4=Quad2D4Node_Stress(E,NU,0.5,0.5,0,0.5,0,0,0.5,0,u4,ID)总体刚度矩阵：各节点位移：各单元应力：三．结果各个节点位移：u1=1.5749，v1=-4.5116，u2=0.5858，v2=-4.2489，u3=-0.4401，v3=-4.1495，u4=1.1458，v4=-3.3911，u5=0.7035，v5=-2.9251，u6=-0.4105,v6=-3.0964，u7=0，v7=-3.0486，u8=0.6532，v8=-1.9914，u9=0，v9=0。单元应力分布：单元1的应力分量为Ϭx=0.1379Pa,Ϭy=-0.6942Pa,τxy=-0.4052Pa;单元2的应力分量为Ϭx=-0.1379Pa,Ϭy=0.0375Pa,τxy=-0.0948pa;单元3的应力分量为Ϭx=0.8696Pa,Ϭy=-1.3058Pa,τxy=-0.5948pa;单元4的应力分量为Ϭx=-0.8696Pa,Ϭy=-2.0375Pa,τxy=-0.9052pa。