有限元作业：三角形单元求解.doc

有限元作业年 级2015级学 院机电工程学院专业名称班级学号学生姓名2016年05月 如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比V=1/6,厚度为I=1。按平面应力问题计算，运用有限元方法，分别采用三角形及四边形单元求解，求节点位移及单元应力（要求三角形单元数量不少于4个，四边形单元不少于2个）图(一)图（二）三角形单元求解图（三）四边形单元求解（1） 如图划分三角形单元，工分成四个分别为（2） 如图分别进行编号1、2、3、4、5、6，并建立坐标系（3） 编程进行求解，得出结果，其中假设力P=2000N调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k1 = 1.0e+06 * 7.2857 -3.0000 -2.1429 0.8571 -5.1429 2.1429 -3.0000 7.2857 2.1429 -5.1429 0.8571 -2.1429 -2.1429 2.1429 2.1429 0 0 -2.1429 0.8571 -5.1429 0 5.1429 -0.8571 0 -5.1429 0.8571 0 -0.8571 5.1429 0 2.1429 -2.1429 -2.1429 0 0 2.1429k2 = 1.0e+06 * 5.1429 0 -5.1429 0.8571 0 -0.8571 0 2.1429 2.1429 -2.1429 -2.1429 0 -5.1429 2.1429 7.2857 -3.0000 -2.1429 0.8571 0.8571 -2.1429 -3.0000 7.2857 2.1429 -5.1429 0 -2.1429 -2.1429 2.1429 2.1429 0 -0.8571 0 0.8571 -5.1429 0 5.1429k3 = 1.0e+06 * 2.1429 0 -2.1429 -2.1429 0 2.1429 0 5.1429 -0.8571 -5.1429 0.8571 0 -2.1429 -0.8571 7.2857 3.0000 -5.1429 -2.1429 -2.1429 -5.1429 3.0000 7.2857 -0.8571 -2.1429 0 0.8571 -5.1429 -0.8571 5.1429 0 2.1429 0 -2.1429 -2.1429 0 2.1429k4 = 1.0e+06 * 2.1429 0 -2.1429 -2.1429 0 2.1429 0 5.1429 -0.8571 -5.1429 0.8571 0 -2.1429 -0.8571 7.2857 3.0000 -5.1429 -2.1429 -2.1429 -5.1429 3.0000 7.2857 -0.8571 -2.1429 0 0.8571 -5.1429 -0.8571 5.1429 0 2.1429 0 -2.1429 -2.1429 0 2.1429调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵求出的节点位移U = 0 0 0 0 -0.0004 0.0008 0.0005 0.0010 0.0007 0.0023 -0.0007 0.0026调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,SxyS1 = 1.0e+03 * -4.4086 -0.7348 3.5914S2 = 1.0e+03 * 4.4086 -0.6405 0.4086S3 = 1.0e+03 * 1.8907 -1.0601 2.1093S4 = 1.0e+03 * -1.8907 2.10931.8907二、（1）如图划分四边形单元，工分成四个分别为（2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N调用 Quad2D4Node_Stiffness函数，求出单元刚度矩阵 调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵求出节点位移U = 0 0 0 0 0.0012 0.0017 -0.0012 0.0017 0.0016 0.0049 -0.0017 0.0052调用Quad2D4Node_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量S1 = 1.0e+03 * 0.0000 -0.2478 2.0000S2 = 1.0e+07 * 0.6856 4.1135 -1.7137程序附录一、1、三角形单元总程序：E=1e7;NU=1/6;t=1;ID=1;%调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k1=Triangle2D3Node_Stiffness(E,NU,t,0,1,0,0,1,1,ID)k2=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,1,1,ID)k3=Triangle2D3Node_Stiffness(E,NU,t,1,1,1,0,2,0,ID)k4=Triangle2D3Node_Stiffness(E,NU,t,2,0,2,1,1,1,ID)%调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵KK = zeros(12,12);KK=Triangle2D3Node_Assembly(KK,k1,1,2,3);KK=Triangle2D3Node_Assembly(KK,k2,2,4,3);KK=Triangle2D3Node_Assembly(KK,k3,3,4,5);KK=Triangle2D3Node_Assembly(KK,k4,5,6,3)% 边界条件的处理及刚度方程求解k=KK(5:12,5:12) p=0;0;0;0;0;0;0;2000u=kp%支反力的计算U=0;0;0;0;u %为节点位移P=KK*U%调用Triangle2D3Node_Strain函数，求出应变SN1、SN2、SN3中求出的分别为SNx,SNy,SNxyu1=U(1);U(2);U(3);U(4);U(5);U(6);u2=U(3);U(4);U(7);U(8);U(5);U(6);u3=U(5);U(6);U(7);U(8);U(9);U(10);u4=U(9);U(10);U(11);U(12);U(5);U(6); SN1=Triangle2D3Node_Strain(0,1,0,0,1,1,u1) SN2=Triangle2D3Node_Strain(0,0,1,0,1,1,u2) SN3=Triangle2D3Node_Strain(1,1,1,0,2,0,u3) SN4=Triangle2D3Node_Strain(2,0,2,1,1,1,u4) %调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,Sxy u1=U(1);U(2);U(3);U(4);U(5);U(6);u2=U(3);U(4);U(7);U(8);U(5);U(6);u3=U(5);U(6);U(7);U(8);U(9);U(10);u4=U(9);U(10);U(11);U(12);U(5);U(6); S1=Triangle2D3Node_Stress(E,NU,0,1,0,0,1,1,u1,ID) S2=Triangle2D3Node_Stress(E,NU,0,0,1,0,1,1,u2,ID) S3=Triangle2D3Node_Stress(E,NU,1,1,1,0,2,0,u3,ID) S4=Triangle2D3Node_Stress(E,NU,2,0,2,1,1,1,u4,ID)2、求刚度矩阵程序function k=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)%该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU，厚度t%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输出单元刚度矩阵k(6X6)%-A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj)/2;betai = yj-ym;betaj = ym-yi;betam = yi-yj;gammai = xm-xj;gammaj = xi-xm;gammam = xj-xi;B = betai 0 betaj 0 betam 0 ; 0 gammai 0 gammaj 0 gammam ; gammai betai gammaj betaj gammam betam/(2*A);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;endk= t*A*B*D*B;3、 求整体刚度矩阵function z = Triangle2D3Node_Assembly(KK,k,i,j,m)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k%输入单元的节点编号I、j、m%输出整体刚度矩阵KK %-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;for n1=1:6 for n2=1:6 KK(DOF(n1),DOF(n2)= KK(DOF(n1),DOF(n2)+k(n1,n2); endendz=KK;4、求应变程序function strain=Triangle2D3Node_Strain(xi,yi,xj,yj,xm,ym,u)%该函数计算单元的应变%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入单元的位移列阵u(6X1)%输出单元的应力strain(3X1)，由于它为常应变单元，则单元的应变分量为SNx,SNy,SNz%-A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj)/2;betai = yj-ym;betaj = ym-yi;betam = yi-yj;gammai = xm-xj;gammaj = xi-xm;gammam = xj-xi;B = betai 0 betaj 0 betam 0 ; 0 gammai 0 gammaj 0 gammam ; gammai betai gammaj betaj gammam betam/(2*A);strain = B*u;5、求应力程序function stress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%该函数计算单元的应力%输入弹性模量E，泊松比NU，厚度t%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)，单元的位移列阵u(6X1)%输出单元的应力stress(3X1)，由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy%-A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj)/2;betai = yj-ym;betaj = ym-yi;betam = yi-yj;gammai = xm-xj;gammaj = xi-xm;gammam = xj-xi;B = betai 0 betaj 0 betam 0 ; 0 gammai 0 gammaj 0 gammam ; gammai betai gammaj betaj gammam betam/(2*A);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;endstress = D*B*u;二、1、四边形单元总程序：E=1e7;NU=1/6;h=1;ID=1;%调用 Quad2D4Node_Stiffness函数，求出单元刚度矩阵 k1= Quad2D4Node_Stiffness(E,NU,h,0,1,0,0,1,0,1,1,ID) k2= Quad2D4Node_Stiffness(E,NU,h,1,0,2,0,2,1,1,1,ID) %调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵 KK=zeros(12,12); KK= Quad2D4Node_Assembly(KK,k1,1,2,3,4); KK= Quad2D4Node_Assembly(KK,k2,3,5,6,4) % 边界条件的处理及刚度方程求解 k=KK(5:12,5:12) p=0;0;0;0;0;0;0;2000u=kp%支反力的计算U=0;0;0;0;u %为节点位移P=KK*U%调用Quad2D4Node_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量u1=U(1);U(2);U(3);U(4);U(5);U(6);U(7);U(8);u2=U(5);U(6);U(9);U(10);U(11);U(12);U(7);(8);S1= Quad2D4Node_Stress(E,NU,0,1,0,0,1,0,1,1,u1,ID)S2= Quad2D4Node_Stress(E,NU,1,0,2,0,2,1,1,1,u2,ID)2、求刚度矩阵程序function k= Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)%该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU，厚度h%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输出单元刚度矩阵k(8X8) %-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 xj xm xp*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);3、求总体刚度矩阵程序function z = Quad2D4Node_Assembly(KK,k,i,j,m,p)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k，单元的节点编号i、j、m、p%输出整体刚度矩阵KK%-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;4、求应力程序function stress= Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)%该函数计算单元的应力%输入弹性模量E，泊松比NU，厚度h，%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp,%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输入单元的位移列阵u(8X1)%输出单元的应力stress(3X1)%由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy%-syms s t;a = (yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s)/4;b = (yi*(t-1)+yj*(1-t)+y