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 = (yi*(t-1)+y

2、j*(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)/4 a*(1-t)/4-

3、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-s s+t -t-1 0;

4、 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); endfunction z = Qu

5、ad2D4Node_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,u,ID) syms s

6、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-t)/4-b*(-1-s)

7、/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; Jfirst = 0

8、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(str1, s,t, 0,0);

9、 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)

10、;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);U(11);U(12)

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

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