有限元课程设计
一. 问题描述
如图所示的平面矩形结构,设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)+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)/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-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]; end
BD = J*transpose(B)*D*B; r = int(int(BD, t, -1, 1), s, -1, 1); z = h*r;
k = double(z); end
function 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); end end z=KK; end
function stress= Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,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-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];
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]; end
str1 = 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\p
U=[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,τ单元2的应力分量为Ϭx=-0.1379Pa,Ϭy=0.0375Pa,τ单元3的应力分量为Ϭx=0.8696Pa,Ϭy=-1.3058Pa,τ单元4的应力分量为Ϭx=-0.8696Pa,Ϭy=-2.0375Pa,τ
xy
=-0.4052Pa; =-0.0948pa; =-0.5948pa;
xy
xy
xy
= -0.9052pa。