有限元作业三角形单元求解文档格式.docx

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、 2.1429 -2.1429 -2.1429 0 0 2.1429k2 = 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 = 2.1429 0 -2.1429 -2.1429 0 2.14

3、29 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 0k4 =调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵求出的节点位移U = 0 -0.0004 0.0008 0.0005 0.0010 0.0007 0.0023 -0.0007 0.0026调用Triangle2D3Node_Stress函数，

4、求出应力,S1、S2、S3、中求出的分别为Sx,Sy,SxyS1 = 1.0e+03 * -4.4086 -0.7348 3.5914S2 = 4.4086 -0.6405 0.4086S3 = 1.8907 -1.0601 2.1093S4 = -1.89071.8907二、（1）如图划分四边形单元，工分成四个分别为调用 Quad2D4Node_Stiffness函数，求出单元刚度矩阵调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵求出节点位移 0.0012 0.0017 -0.0012 0.0016 0.0049 -0.0017 0.0052调用Quad2D4Node

5、_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量 0.0000 -0.2478 2.0000 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_Stif

6、fness（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）%

7、边界条件的处理及刚度方程求解k=KK（5:12,5:12） p=0;0;2000u=kp%支反力的计算U=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（7）;U（8）;u3=U（5）;U（6）;U（9）;U（10）;u4=U（9）;U（10）;U（11）;U（12）; SN1=Triangle2D3Node_Strain（0,1,0,0,1,1,u1） SN2=Triangle2D3Node_Stra

8、in（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）; 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

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

10、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

11、 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: KK（DOF（n1）,DOF（n2）= KK（DOF（n1）,DOF（n2）+k（n1,n2）; endz=KK;4、求应变程序function strain=Triangle2D3Node_Strain（xi,yi,xj,yj,xm,ym,u）%该函数计算单元的应变%输入单元的位移列阵u（6X1）%输出单元的应力strain（3X1），由于它为常应变单元，则单元的应变分量为SNx,SNy,SNz%-