有限元课程设计

1、精选优质文档-倾情为你奉上一 问题描述如图所示的平面矩形结构，设E=1，NU=0.25，h=1，考虑以下约束和外载：位移边界条件BC（u）：UA=0，VA=0，UD=0，力边界条件BC（p）：在CD边上有均布载荷q=1，建模情形：使用四个四节点矩形单元，试在该建模情形下，求各节点的位移以及各个单元的应力分布。二 Matlab程序(1).函数定义：function k= Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID) syms s t; a = (yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s)/4; b

2、= (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)/4 0 ; 0 c*(s-1)/4-d*(t-1)/4 ; c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4; B2 = a*(1-t)/4-b*(-1-s)/4 0 ; 0 c*(-1-s)/4-d*(1-t)/4 ; c*(-1-s)/4-d*(1-t)

3、/4 a*(1-t)/4-b*(-1-s)/4; B3 = a*(t+1)/4-b*(s+1)/4 0 ; 0 c*(s+1)/4-d*(t+1)/4 ; c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4; B4 = a*(-1-t)/4-b*(1-s)/4 0 ; 0 c*(1-s)/4-d*(-1-t)/4 ; c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4; Bfirst = B1 B2 B3 B4; Jfirst = 0 1-t t-s s-1 ; t-1 0 s+1 -s-t ; s-t -s-1 0 t+1 ; 1-

4、s s+t -t-1 0; J = xi xjxmxp*Jfirst*yi ;yj ; ym ; yp/8; B = Bfirst/J; if ID = 1 D = (E/(1-NU*NU)*1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2; elseif ID = 2 D = (E/(1+NU)/(1-2*NU)*1-NU NU 0 ; NU 1-NU 0 ; 0 0 (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); endfu

5、nction z = 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; for n1=1:8 for n2=1:8 KK(DOF(n1),DOF(n2)= KK(DOF(n1),DOF(n2)+k(n1,n2); endendz=KK;endfunction stress= Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,

6、u,ID) syms s t; 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)/4 0 ; 0 c*(s-1)/4-d*(t-1)/4 ; c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4; B2 = a*(1-

7、t)/4-b*(-1-s)/4 0 ; 0 c*(-1-s)/4-d*(1-t)/4 ; c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4; B3 = a*(t+1)/4-b*(s+1)/4 0 ; 0 c*(s+1)/4-d*(t+1)/4 ; c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4; B4 = a*(-1-t)/4-b*(1-s)/4 0 ; 0 c*(1-s)/4-d*(-1-t)/4 ; c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4; Bfirst = B1 B2 B3 B4

8、; Jfirst = 0 1-t t-s s-1 ; t-1 0 s+1 -s-t ; s-t -s-1 0 t+1 ; 1-s s+t -t-1 0; J = xi xjxmxp*Jfirst*yi ;yj ; ym ; yp/8; B = Bfirst/J; if ID = 1 D = (E/(1-NU*NU)*1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2; elseif ID = 2 D = (E/(1+NU)/(1-2*NU)*1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2; endstr1 = D*B*u; str2 = subs(str

9、1, 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,

10、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=kpU=u(1:12);0;u(13:15);0;0;u1=U(1);U(2

11、);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,

12、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。